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ABSTRACT 

In  this  paper  we  present  a  class  of  ENO  schemes  for  the  numerical  solution  of  multidimen¬ 
sional  hyperbolic  systems  of  conservation  laws  in  structured  and  unstructured  grids.  This 
is  a  class  of  shock-capturing  schemes  which  are  designed  to  compute  cell-averages  to  high- 
order  of  accuracy.  The  ENO  scheme  is  composed  of  a  piecewise-polynomial  reconstruction 
of  the  solution  from  its  given  cell-averages,  approximate  evolution  of  the  resulting  initial- 
value  problem,  and  averaging  of  this  approximate  solution  over  each  cell.  The  reconstruction 
algorithm  is  based  on  an  adaptive  selection  of  stencil  for  each  cell  so  as  to  avoid  spurious 
oscillations  near  discontinuities  while  achieving  high  order  of  accuracy  away  from  them. 


’Research  was  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Contract 
No.  NAS1-18605  while  the  author  was  in  residence  at  the  Institute  for  Computer  Applications  in  Science 
and  Engineering  (ICASE),  NASA  Langley  Research  Center.  Hampton,  VA  23665.  Research  also  supported 
by  ONR  grant  N00014-86-K-0691,  DARPA  grant  in  the  ACMP  program  and  NSF  grant  DMS  88-11863. 
t  Research  supported  by  Rockwell  International  Science  Center’s  IRckD  funds. 
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1.  Introduction 


In  this  paper  we  generalize  the  ENO  schemes  of  [8],  [9]  and  [3].  to  several  space 
dimensions  with  structured  and  unstructured  grids;  it  should  be  read  in  conjunction 
with  these  earlier  articles. 

The  ENO  schemes  are  of  the  form 

un+1  =  A  ■  E(t)  ■  R(-,  vn )  =  Eh-vn  ,  (1.1) 

where  vn  are  cell-averages  of  u(x,t„),  the  solution  at  time  tn.  jF2(x;  vn)  is  a  reconstruc¬ 
tion  procedure  which  produces  a  high-order  accurate  global  approximation  to  u(x,tn) 
from  its  given  cell-averages  vn ;  in  this  paper  we  consider  R  to  be  a  polynomial  function 
in  each  cell.  E(t)  is  the  evolution  operator  of  the  PDE  which  includes  the  influence  of 
boundary  conditions.  A  is  the  cell-averaging  operator.  When  we  consider  E(t)  which 
corresponds  to  a  PDE  with  divergence  free  form,  the  scheme  Eh  (1.1)  is  automatically 
in  conservation  form  no  matter  what  is  the  particular  shape  of  the  cells. 

The  scheme  (1.1)  is  linked  to  a  grid  in  a  very  loose  way.  Rewriting  (1.1)  as 

Eh  =  R- A- E  =  Ph  ■  E  (1.2) 

we  can  view  the  scheme  as  a  composition  of  exact  evolution  E  with  a  projection  Ph  = 
R  •  A ,  the  role  of  which  is  to  project  the  solution  into  a  finite-dimensional  space  of 
functions.  The  averaging  A  and  the  reconstruction  R  may  even  use  a  different  set  of 
cells.  This  observation  is  particularly  useful  for  purposes  of  various  grid  manipulations 
like  component  grids,  multigrid  calculations,  and  time-dependent  adaptive  grids. 

The  question  of  stability  and  convergence  of  numerical  schemes  is  related  to  the 
boundedness  or  possible  growth  of  spurious  oscillations  in  the  computed  solution.  The 
largest  source  of  spurious  oscillations  in  the  numerical  solution  is  a  Gibbs-like  phe¬ 
nomenon  associated  with  interpolation  through  a  discontinuity;  these  are  0(1)  oscilla¬ 
tions  with  respect  to  refinement. 
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The  particular  form  of  the  scheme  (1.1)  leaves  the  question  of  control  over  spurious 
oscillations  to  the  design  of  the  reconstruction  i?,  which  is  a  problem  in  the  approxima¬ 
tion  of  functions.  The  ENO  schemes  attempt  to  avoid  growth  of  spurious  oscillations 
by  an  adaptive-stencil  approach,  in  which  each  cell  is  assigned  its  own  stencil  of  cells  for 
purposes  of  reconstruction.  For  each  cell  we  select  an  interpolating  stencil  in  which  the 
solution  is  smoothest  in  some  sense.  Thus  cells  near  a  discontinuity  are  assigned  sten¬ 
cils  from  the  smooth  part  of  the  solution  and  a  Gibbs-like  phenomenon  is  so  avoided. 
The  term  essentially  non-oscillatory  is  used  because  spurious  oscillations  on  the  scale 
of  the  interpolation  error  in  the  smooth  part  of  the  solution  are  not  ruled  out.  This 
adaptive-stencil  strategy  seems  to  ensure  the  stability  and  convergence  of  the  scheme 
0-1). 

The  question  of  accumulation  of  error,  i.e.,  the  relation  between  the  local  trunca¬ 
tion  error  and  the  actual  error  at  the  end  of  the  calculation  is  related  to  the  nature  of 
the  stability  of  the  scheme.  In  [7]  we  have  examined  a  problem  in  which  the  selection 
of  smoothest  stencil  for  the  reconstruction  leads  to  a  scheme  (1.1)  which  is  linearly 
unstable  in  the  whole  interval.  We  have  found  that  once  the  high-order  derivatives 
have  begun  to  oscillate  and  grow,  the  selection  of  stencil  became  erratic  and  this  has 
stabilized  the  computation,  so  the  calculation  was  convergent  but  with  a  reduced  order 
of  accuracy.  Later  Meiburg  [13]  has  shown  that  similar  loss  of  accuracy  can  occur  near 
a  point  at  which  all  the  derivatives  up  to  a  certain  order  vanish.  Shu  [14]  has  demon¬ 
strated  that  this  unnecessary  loss  of  accuracy  can  be  avoided  by  biasing  the  selection 
of  the  stencil  toward  a  central  one  in  the  smooth  part  of  the  solution  —  we  adopt  this 
strategy  here  as  well. 

In  Section  2  we  introduce  notation  which  enables  us  to  describe  the  scheme  (1.1) 
in  the  most  general  case.  In  Section  3  we  describe  the  process  of  reconstruction  in  a 
given  cell.  After  these  preliminaries  we  begin  to  tackle  the  two  really  important  issues 
of  designing  a  procedure  for  selecting  a  stencil  in  general  geometries  and  the  practical 
aspect  of  an  efficient  implementation. 
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In  Section  4  we  go  through  several  levels  of  approximation  for  the  numerical  flux, 
to  finally  arrive  at  a  simple  expression  which  is  easy  to  use  and  is  yet  adequate. 

In  Section  5  we  describe  two  efficient  implementations  of  the  resulting  scheme. 

In  Section  6  we  examine  the  schemes  (1.1)  with  a  fixed  central  stencil  which  is 
to  be  used  with  some  form  of  added  numerical  dissipation.  In  this  context  we  outline 
a  new  technique  to  obtain  an  ENO  reconstruction  within  the  fixed  central  stencil  by 
hybridizing  the  high-order  reconstruction  with  the  first-order  one.  In  Appendix  1  we 
present  more  details  and  analysis.  We  bring  this  preview  of  a  future  paper  here  in  older 
to  give  a  general  picture. 

In  Section  7  we  present  the  adaptive  algorithm  for  the  selection  of  an  appropriate 
stencil  of  cells  in  the  most  general  case  and  describe  its  efficient  implementation  as  an 
ordered  Gauss  elimination  with  adaptive  row-pivoting.  At  the  end  of  this  section  we 
show  how  to  use  this  procedure  to  automatically  select  fixed  stencils  which  are  either 
central  or  directionally  biased.  We  also  show  how  to  bias  the  ENO  stencil  toward  a 
central  one. 

In  Section  8  we  discuss  special  reconstruction  techniques  for  solutions  of  hyperbolic 
systems  of  conservation  laws.  These  special  techniques  may  be  needed  in  order  to  handle 
particularly  strong  interaction  of  discontinuities  in  the  computed  solution. 

In  Section  9  we  describe  the  application  of  ENO  schemes  for  rectangular  _,nds  and 
show  the  great  simplification  that  occurs  in  this  case.  For  rectangular  and  smoothly 
varying  grids  we  use  reconstruction  via  deconvolution  in  order  to  obtain  efficient  schemes 
that  use  a  tensor-product  of  one- dimensional  stencils.  The  third-or  jer  accurate  scheme 
turns  out  to  be  particularly  simple,  and  we  feel  that  it  is  of  immediate  practical  impor¬ 
tance  as  the  next  generation  to  the  second-order  accurate  T\  D  schemes.  In  Appendix  2 
we  describe  the  implementation  of  the  third-order  scheme  to  the  solution  of  the  Euler 
equations  of  gasdynamics. 
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In  Section  10  we  consider  the  question  of  time-integration  and  compare  the  relative 
merits  of  a  method  of  lines  approach  to  a  local  Cauchy-Kowalewski  procedure  of  a  single 
step.  The  method  of  lines  is  easier  to  program  but  more  expensive  to  use.  Thus  we 
advocate  the  use  of  the  one-step  procedure  for  production  codes  whenever  feasible. 

In  Appendix  3  we  outline  a  new  scheme  which  uses  alternating  dual  grids.  The 
two  sets  of  grids  correspond  to  centroids  and  vertices.  The  new  scheme  alternatingly 
computes  point-values  in  one  set  of  cells  and  cell-averages  in  the  dual  one,  at  no  extra 
computational  cost.  The  use  of  these  two  sets  of  values  enables  us  to  obtain  a  more 
compact  reconstruction.  We  present  this  “2  for  the  price  of  1”  scheme  because  we  feel 
that  this  is  the  next  step  in  the  development  of  high-order  ENO  schemes. 

In  writing  this  paper  we  have  attempted  to  give  as  broad  a  picture  as  possible  on 
the  development  of  ENO  schemes  as  we  see  it.  The  number  of  schemes  which  can  be 
formed  by  different  numerical  fluxes,  various  reconstruction  techniques,  and  methods 
of  time-integration  is  huge.  Consequently  our  numerical  experimentation  amounts  to  a 
mere  sampling.  The  use  of  analysis  for  the  design  of  these  highly  nonlinear  schemes  is 
still  very  limited  and  one  has  to  rely  heavily  on  numerical  experimentation.  We  hope 
that  this  article  will  encourage  others  to  experiment  with  these  schemes. 

We  would  like  to  point  out  some  related  work  that  we  know  of.  The  scheme  (1.1) 
was  originated  by  Godunov  to  design  his  first-order  scheme  [6],  and  was  subsequently 
extended  to  second-order  accuracy  by  van  Leer  [17],  and  improved  upon  by  Colella  and 
Woodward  [4], 

The  case  r  =  2  of  the  ENO  schemes  also  corresponds  to  second-order  accurate 
TVD  schemes;  thus  all  finite- volume  TVD  schemes  are  related  to  it.  In  this  context  we 
would  like  to  refer  the  reader  to  a  recent  work  by  Durlofsky,  Osher,  and  Engquist  [5] 
on  second-order  TVD  schemes  in  a  triangular  grid. 

We  would  also  like  to  point  out  the  work  of  Casper  [2]  on  fourth-order  accurate 
ENO  implementation  for  rectangular  and  smoothly  varying  grids  and  the  work  of  Shu 
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and  Osher  [15].  The  latter  scheme  is  an  efficient  ENO  scheme  for  point-values  on  rect¬ 
angular  and  smoothly  varying  grids,  and  thus  is  not  of  the  form  (1.1).  The  conservation 
form  of  this  scheme  is  obtained  by  a  clever  trick  in  which  the  numerical  flux  is  treated 
as  a  primitive  of  some  other  function. 

We  would  also  like  to  refer  the  reader  to  a  recent  paper  by  Barth  and  Fredrickson 
[1]  who  have  implemented  the  scheme  (1.1)  for  unstructured  triangular  grids  using 
large  fixed  central  stencils  with  least  squares  reconstruction.  Their  numerical  results 
demonstrate  that  a  high  order  of  accuracy  can  be  achieved  for  smooth  solutions  even 
on  highly  irregular  grids. 
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2.  The  Numerical  Scheme 


In  this  paper  we  consider  the  Initial  Boundary  Value  Problem  (IBVP)  for  a  hy¬ 
perbolic  system  of  conservation  laws  in  S-space  dimensions: 

ut  +  div  f(u)  =  0  ,  x  G  D  C  TZS  ,  t  >  0  (2.1a) 

u(x,  0)  =  uq(x)  ,  x  €  T>  (2.1b) 


with  given  boundary  conditions  on  dV.  the  boundary  of  V.  The  computational  domain 
V  is  divided  into  cells  Cj 

X>=U  C>  •  C.flC^O  .  (2.2a) 

We  assume  that  dCj ,  the  boundary  of  the  cell  Cj ,  is  piecewise-smooth,  i.e., 


8Cj  =  |J  dC)  (2.2b) 

k 

where  dCj  is  smooth.  Usually  OCj  is  linear,  yet  our  formulations  allow  for  nonzero 
curvature.  We  also  assume  that  there  is  a  refinement  parameter  h  such  that  the  largest 
sphere  contained  in  the  cells  is  of  radius  0(h) ,  and  that  the  ratio  between  the  largest 
sphere  to  the  smallest  one  in  the  computational  domain  remains  bounded  under  refine¬ 
ment.  We  denote 


(2.2c) 


and  by  Cj,  the  centroid  of  the  cell  C} 


c,  = 


(2. 2d) 


Let  Uj  denote  the  cell-average  of  the  function  u(x )  over  Cj 

Uj  =  —7-7  [  u(x)dx  =  A(Cj)u  (2.3) 

Il  jI  Jr, 

and  denote  by  A(C} )  the  cell- averaging  operator.  Given  cell-averages  u  =  { IZj }  of  u(x) 
in  P,  we  denote  by  R(x;u  )  a  reconstruction  of  u  from  u  which  satisfies 

R(x;u  )  =  u(.r)  4-  O  (hr )  wherever  u  is  smooth  (2.4a) 

A  [Cj  )/?(•;  u  )  =  Uj  (conservation)  .  (2.4b) 
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Typically  R  is  a  piecewise- polynomial  function  of  degree  r  —  1  which  is  discontin¬ 
uous  across  the  cell-houndary  dCj.  Let  Rj(x-,vn)  denote  the  polynomial  which  defines 
R  in  Cj,  i.e., 

R{x-v11)  =  Rj(x-vn)  for  x  e  Cj  (2.4c) 

Let  E(t )  denote  the  evolution  operator  of  the  IBVP  (2.1) 

n(-,t)^E(t)u0  (2‘5) 

Note  that  E(t)  includes  the  boundary  conditions  on  dT>. 

We  turn  now  to  describe  the  numerical  scheme  which  is  an  explicit  method  for  the 
approximation  of  the  cell-averages  of  u(x,  t ) 

v”  «  A  (Cj)  «(•,<„)  (2.6) 

We  initialize  the  computation  by  setting 

v)  =  A{Cj)uQ  (2.7a) 

where  u 0  is  the  initial  value  (2.1b).  Given  vn  =  {u"  }  we  compute  v"+1  by 

u"+1  =A(C})E(r)R(-,vn)  (2.7b) 

Thus  we  first  get  a  piecewise-polynomial  approximation  i?(x;  vn)  to  the  solution  u(x,tn). 
Then  we  apply  E(t)  to  R(x;un),  i.e..  we  get  a  solution  in  the  small  (small  r)  to  the 

IBVP 


ic(  +  div  /( tv)  =  0  ,  x  £  V  ,  0<f<r  (2.8a) 

u>(x,0)  =  R(x;  vn)  ,  (2.8b) 

with  the  given  boundary  conditions  on  &D. 

Integrating  the  PDE  (2.8a)  over  C,  x  [0,r]  and  using  the  divergence  theorem  we 
get 

m  - »;)  +  ['  i  [/(u.(t,o)  •  m  ds  =  o  (2. Sc) 

Jo  JdCj 
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Here  N  is  the  outward  normal  to  dCj  and  the  integration  is  over  the  boundary  of  the 
cell.  Hence  the  application  of  A(Cj)  to  the  solution  of  (2.8a),  (2.8b)  at  time  r  is  given 
by 

»;+I  1^7  f  {  [f(m  •  R(-;v"))  ■  IVJdS  .  (2.9) 

IWI  Jo  Jac j 

The  above  expression  is  an  explicit  scheme  in  conservation  form.  We  point  out  that 
the  term  vj  in  (2.9)  is  obtained  from 

A(Cj)w{x,  0)  =  A(Cj)R(x-vn)  =  vj  .  (2.10) 

Thus  property  (2.4b)  of  the  reconstruction  is  essential  in  order  to  get  conservation  form 
from  (2.7b). 


Since  A(Cj)  is  a  positive  operator  and  £(r)  is  the  exact  evolution  operator,  the  con¬ 
trol  over  possible  growth  of  oscillations  in  the  numerical  solution  is  applied  through  the 
reconstruction  R(x-,  vn).  In  Section  7,  we  shall  describe  an  Essentially  Non-Oscillatory 
(ENO)  reconstruction  technique  which  is  designed  in  order  to  achieve  this  goal. 


In  this  paper  we  concentrate  on  the  semi-discrete  formulation  of  (2.9),  which  is 
obtained  by  dividing  (2.8c)  by  r  and  letting  r  — >  0.  In  all  the  cells  which  do  not  have 
a  common  side  with  the  boundary  dD  we  get 


(2.11a) 


In  (2.11)  we  have  introduced  some  new  notations  and  conventions: 


In  =  /  •  N 


(2.11b) 


/$  (u-i,u2)  =  /n  (^r(0;u1,u2)) 


(2.11c) 


where  W  (£/t;  Hi ,  u2)  is  the  self-similar  solution  of  the  one- dimensional  Riemann  prob¬ 


lem 


wt  -f  ^7/n{w)  =  0 

ot, 


u>(£,  o) 


f  «1  (.  <  o 
\  £  >  o 


(2.12a) 

(2.12b) 
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Observe  that  £  is  scalar.  Rj  in  the  term  is  the  polynomial  (2.4c)  which 

is  evaluated  at  the  boundary  dCy  R »  is  the  reconstruction  in  the  cell  which  is  in  the 
exterior  of  Cj  and  has  dCj  in  common.  Thus  *  is  a  symbolic  notation  for  the  index  of 
such  a  cell.  Another  innovation  is  the  breakup  of  the  integral  over  dC}  into  its  smooth 
pieces  dC '*  in  (2.2b).  (The  superscript  R  in  fj$  stands  for  “Riemann”.) 


When  dCj  is  an  element  of  the  boundary  dV  we  make  the  following  substitution 
in  (2.11) 

/  ffj  (Rj,  f?*)  ds  — >  [  fN(E(0+)Rj)ds  (2.1T) 

iac*  Jacf 

which  is  computed  at  dT>  and  includes  the  boundary  conditions. 
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3.  Reconstruction 


Given  cell-averages  u  =  {ttj}  of  u(x),  we  describe  in  this  section  a  reconstruction 
technique  I ?(x:n  )  which  satisfies  properties  (2.4) 

R(.r\u  )  =  u(x)  +  0  (h' )  wherever  u(.r)  is  smooth  (accuracy)  (3.1a) 
A(Cj)R{-\u  )=  Tij  (conservation)  .  (3.1b) 


In  this  paper  we  consider  a  piecewise-polynomial  reconstruction,  which  is  defined 
by  a  polynomial  of  degree  ( r  —  1)  in  each  of  the  cells.  As  in  (2  4c)  we  denote  the 
polynomial  in  the  cell  C,  by  R,(x-,u  )  and  express  it  as  a  Taylor  expansion  around  the 
centroid  c^: 

r_1  1 

R{x:u)  =  Ri(x-,n)  =  ^2—^2(x-cl)tDe  ,  x  €  C,  (3.2) 

k= 0  k'  |l|=& 

Here  we  use  a  multi-index  notation  and  convention 


£  —  (G i • '  •  ■  O  ,  \t\-£i+ (2 +  ■••  +  £, 

ye  =  (yi)/|  •(!&)**  •••(y.)/* 

The  summation  convention  in  stands  for 

k  k  k 

E  =  EE-E 

\(\=zk  ti=0  t2=o  t,= 0 

b  +  (2  1 - b  =  1' 


(3.3a) 

(3.3b) 


(3.3c) 


so  that  terms  D(  corresponding  to  mixed  derivatives  appear  in  the  summation  ex¬ 
actly  the  number  of  times  they  should.  We  also  use  the  multi-index  convention  for 
differentiation 


df  _  0W 

dx[l  dx(27  •  •  ■  dx(s’ 


(3.3d) 


Comparing  (3.2)  with  the  Taylor  expansion  of  u(x)  around  the  centroid  C, 


deu 


,  k\ 

k~  0  \(\z=k 


dx( 


Ic, 


+  0(hr) 


(3.4) 
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we  see  that  the  accuracy  requirement  (3.1a)  amounts  to 


D,  = 


(a)  +  o(hr-'e ,  0  <  |f|  <  r  —  1 


(3.5a) 


We  note  that  (3.5a)  for  |£|  =  0  reads 


Do  =  Ri  (r,;  u  )  =  u  (ci)  +  0(hr)  . 


(3.5b) 


Let  J(i )  be  a  set  of  indices  of  cells  which  includes  i.  We  shall  refer  to  J(i)  as  a 
stencil  of  cells  associated  with  the  cell  i  and  denote  the  number  of  cells  in  it  by  |  J\.  We 
:on«ider  now  the  set  of  linear  equations  for  Dt,  0  <  |f|  <  r  —  1,  wbirh  is  obtained  by 
taking  a  cell-average  of  the  polynomial  Rx(x;u  )  in  (3.2)  over  all  the  cells  j  in  J(i) 


A(Cj)Ri  =  Uj  ,  jeJ(i) 


(3.6a) 


where 


y  y  <*>,*£*  =  «_,•  ,  j  e  j(o 

k- 0  \(\=k 


aht  =  ~A(Cj)(x  -  c,)'  =  - 


dV  . 


To  get  proper  scaling  let  us  also  consider  an  alternative  form  of  (3.6) 

X]  a']fDe  =  u) 

*= 0  |/|=Jt 


whjre 


a'ht  = 

D't  =  h^Di  . 

Next  we  rewrite  this  system  of  linear  equations  in  matrix  form 


(3.6b) 


(3.6c) 


(3.6b') 


Qd  =  u 


(3.7a) 


Here  d  is  the  vector  of  unknowns 


J 

d  =  {d\ ,  •  •  • ,  cbc)  ,  k  —  n(r) 


(3.7b) 
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in  which  D'f  are  arranged  in  groups  of  equal  \i\  with  increasing  order  of  |£|.  Q  is  a  matrix 
with  |  J(i)\  rows  and  n(r)  columns;  let  us  agree  that  the  first  row  always  corresponds  to 
j  =  i.  The  entries  of  Q  are  the  proper  rearrangement  of  the  coefficients  a'-,  in  (3.6b'). 
u  is  the  vector  of  the  given  cell-averages 


=  (»].•  ■  •  ,U|7|)  ,  U]  -  Ui  . 


(3.7c) 


When  we  apply  the  cell-averaging  operator  to  u(x)  in  (3.4)  we  get 


u}  ^  A  (Cj)u(x)  =  +  0(hr) 


k= o  \e\=k 


=  EE'ihf'‘''l0(c’)l+o(',r)  • 

k=0  |/|=Ar  J 


(3.8a) 


Let  us  denote  by  dE  the  vector  (3.7b)  in  which  we  substitute  D'f  by  /i^l^Kc,);  clearly 


QdE  =  u  +  0(hr) 


(3.8b) 


Subtracting  (3.8b)  from  (3.7a)  we  get 


Q{d-dE)  =  0{hr)  . 


(3.8c) 


Let  us  denote  symbolically  by  Q  the  solution  procedure  we  are  going  to  use  in  order 


to  obtain  d  from  a,  i.e., 


d  =  Q  u  . 


(3.9a) 


Q  ||  <  constant  as  h  — ♦  0 


(3.9b) 


then  it  follows  from  (3.8c)  that 


d-dE\\=0(hr)  . 


(3.9c) 


This,  due  to  the  scaling  (3.6b'),  implies  the  accuracy  requirements  (3.5a)  and  (3.1a). 


We  turn  now  to  consider  the  requirement  (3.9b)  from  the  point  of  view  of  the 
stencil  J{i).  Clearly  we  need 

rank  Q  =  «(r)  ,  (3.10a) 

i.e.,  the  number  of  linearly  independent  equations  in  (3.6)  should  be  the  same  as  the 
number  of  terms  in  the  Taylor  expansion  (3.2).  This,  of  course,  implies  that  | «7(?)|,  the 
number  of  cells  in  the  stencil  J(i),  should  be  at  least  /c(r) 

jj(i)|  >  n(r )  .  (3.10b) 

Moreover,  to  assure  that  there  are  enough  linearly  independent  equations  in  the  stencil, 
the  stencil  should  be  large  enough  in  all  spatial  directions,  so  that  all  derivatives  can 
be  properly  approximated  to  the  required  order  of  accuracy  (3.5a). 


When  the  computational  cells  are  defined  by  some  structured  grid,  it  is  possible  to 
predetermine  proper  stencils  with  |  J(«)|  =  «(r).  However,  for  completely  unstructured 
grids  it  seems  wise  to  start  with  very  large  stencils, 

|7(i)|  >>  /c(r)  (3.11a) 

and  to  either  use  a  least-squares  approach,  i.e.,  to  solve  the  k  x  k  system 

QTQd  =  QTu  (3.11b) 


thus  minimizing 


II  Qd~u  ||is 


(3.11c) 


or  to  use  techniques  that  select  «(r)  linearly  independent  equations  out  of  the  many 
available  in  the  stencil.  In  this  latter  category  we  can  use  an  ordered  Gaussian  elimina¬ 
tion  with  row  pivoting  (of  the  type  that  is  described  in  Section  7  for  ENO  reconstruc¬ 
tion)  or  to  get  the  exact  number  of  linearly  independent  equations  needed  by  grouping 
several  of  the  cells  in  J(i)  into  a  single  “super-cell”  and  replace  the  several  equations 
for  the  individual  cells  by  a  single  one  for  the  average  over  the  super-cell.  In  all  such 
techniques  the  process  of  selection  should  be  ordered  so  as  to  give  preference  to  closest 
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neighbors,  i.e.,  the  stencil  should  be  as  centered  around  C,  as  permissible  by  the  various 
constraints. 

We  turn  now  to  examine  the  conservation  property  (3.1b),  and  observe  that  this 
is  exactly  the  equation  for  j  =  i  in  the  system  of  linear  equations  (3.6)  (which  by 
agreement  is  the  first  equation  in  (3.7)).  It  follows  that  if  we  use  the  above  selection 
procedures,  this  property  is  satisfied  automatically.  However,  if  we  use  the  least-squares 
approach  (3.11),  we  have  to  replace  the  first  element  d\  in  the  computed  sohition  by 

n(r) 

d\  =  uj  -  y,  Th  kdk  ,  (3. lid) 

k=\ 

in  other  words,  we  make  sure  that  the  first  equation  is  satisfied  exactly.  In  this  case, 
the  first  equation  may  first  be  eliminated  from  every  other  equation  forming  matrix  Q 
of  rank  n(r)  —  1  and  a  reduced  set  of  equations  may  be  solved  before  applying  (3. lid). 

QrQd  =  Qtu.  (3.11b') 

More  generally,  if  we  compute  approximations  Di  to  the  derivatives  for  \C\  >  2 
which  satisfy  the  accuracy  requirement  (3.5a),  we  define 

r— 1 

Do  =  u.  -  ^  ^2,  •  (3-12) 

k=2\t\=k 

This  ensures  that  the  resulting  reconstruction  (3.2)  satisfies  the  conservation  require¬ 
ment  as  well  as  the  accuracy  requirement  in  (3.1). 

We  point  out  that  the  summation  in  (3.12)  starts  with  k  =  2.  The  reason  for  that 
is  that 

aitt  =  A(C,)  •  (x  -  aY  =  0  for  \i\  =  1  (3.13a) 

due  to  the  fact  that  the  centroid  c,  (2. 2d)  is  defined  by 

c,  =  A{Ci)  x  . 
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It  follows  therefore  that  the  cell-average  is  an  O  ( h 2)  approximation  to  the  point-value 
at  the  centroid 

ut  —  u  ( c, )  =  O  (h2)  .  (3.13b) 

We  remark  that  the  coefficients  in  the  0  (h2)  term  above  are  a,  (  for  |f |  =  2,  which 
depend  on  the  shape  of  the  cell.  Consequently,  numerical  differentiation  of  u ,  (instead 
of  the  point  values  u(c,))  can  give  an  O  (h2)  approximation  to  the  derivatives  of  a, 
only  in  the  case  of  a  smoothly  carving  structured  grid. 

We  end  this  section  with  the  analogous  problem  for  point-values:  Given  point- 
values  of  u(x )  in  the  centroids,  uj  =  u{cj)  f°r  j  €  J(i ),  find  an  rth  order  polynomial 
approximation  /(. r;0)  to  u(x)  in  C,.  Rewriting  this  polynomial  as  a  Taylor  expansion 
around  c, 

r- 1  I 

I(x;u)  =  J2  m  E  (*-Ci)‘Dt  •  (3.14) 

k-  \e\=k 

we  consider  the  system  of  linear  equations  for  Di 

I(cj\u)  =  v.j  ,  j  €  J(i)  (3.15a) 


or  by  the  same  selection  techniques  that  were  mentioned  before. 
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We  observe  that  for  the  reconstruction  (3.Gc)  is  the  cell  average  of  (x  —  c,)#, 
while  aj't  in  (3.15c)  is  the  point-value  of  it.  Consequently, 

Q-><2as/i->0  .  (3.18) 

Hence  for  h  sufficiently  small  we  can  use  the  set  of  cells  which  yield  linearly  independent 
equations  in  Q  for  Q  and  vice  versa. 
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4.  Numerical  Flux 


The  semi-discrete  scheme  (2.11)  can  be  rewritten  in  the  forr 


®vj  _  _  1  dCj  ~fk 

dt  ~  4-  1(7*1  Jh " 


(4.1a) 


where  the  summation  in  k  corresponds  to  the  partition  (2.2b)  of  the  cell  boundary  dCj, 


Here  we  use  the  notation 


7"=wr\L>fti{WS 


~>k  I  f 

j  Jac * 


dS  . 


(4.1b) 


(4.2) 


_ £ 

We  refer  to  the  term  fj  „  in  (4.1b)  as  the  numerical  flux  at  <9(7*;  it  expresses  the 
average  normal  flux  across  dC *  in  the  solution  of  the  IBVP  (2.8).  In  the  following  we 
propose  various  approximations  to  (4.1b)  which  allow  either  a  simplified  or  a  compu¬ 
tationally  more  efficient  algorithm,  or  both.  To  simplify  our  notation,  we  shall  retain 

— k 

the  notation  „  to  the  approximate  as  well  as  the  exact  numerical  flux. 

The  most  straightforward  simplification  is  to  replace  the  integral  in  (4.1b)  by  an 
appropriate  numerical  quadrature, 

/>,.  =  (Rj(* •m\vn),Rt(xm;vn))  (4.3) 

m 

where  xm  are  the  quadrature  points  on  dCj  ,  and  am  are  the  corresponding  quadrature 
coefficients.  The  quadrature  formula  should  be  exact  for  a  polynomial  of  degree  (r  —  1) 
or  more.  Gaussian  quadrature  seems  to  be  particularly  attractive  in  this  context. 

The  next  level  of  simplification  is  to  replace  the  flux  of  the  exact  solution  to  the  Rie- 
mann  problem  (2.11c)  by  a  simpler  approximate  one.  We  observe  that  W  (x/(\ui,u  2), 
the  solution  to  the  Riemann  problem  (2.12),  is  a  Lipschitz  continuous  function  of  uj 
and  U2-  Consequently, 

fs  (ui,u2)  =  ^  [/n(uj)  +  f n{u2 )]  +  0(\u2  -  1/ 1  j )  . 


(4.4a) 


Since 


\Rj{Xm\  vn)  -  i?*(a;m;un)|  =  0{hr) 


(4.4b) 


in  regions  of  smoothness,  we  can  replace  (2.11c)  by 

/jv  («1  >  «2 )  =  \  f//v(ui)  +  fs(u  2)  -  v(ui,u2)(u2  -  u}  )]  (4.5a) 

where  u  is  bounded,  without  adversely  affecting  the  formal  order  of  accuracy  of  the 
scheme.  Again  to  simplify  notation  we  retain  /$  for  the  approximate  fluxes  in  (4.5a). 
The  following  expressions  for  u(u\,u2)  are  suitable: 

v(ui,u2)  =  \As(ui,u2)\  (4.5b) 


where 


An(u,u) 


dfs 

du 


As  (u\,u2)  can  be  taken  as  As  (a),  where  u  is  some  average  of  u  1  and  u2. 


=  |ajv(ui,ti2)|  (4.5b') 

where  as  is  the  maximal  eigenvalue  of  As  (ui,u2). 

v  (ui ,  u2)  =  imix  ja^(v")|  .  (4.5b”) 

The  last  quantity  is  constant  in  space  during  a  time-step,  and  is  computed  anyway  in 
order  to  calculate  the  permissible  time-step  under  a  CFL  restriction. 


We  remark  that  in  the  context  of  first  order  schemes,  (4.5b)  corresponds  to  Roe’s 
scheme,  (4.5b')  corresponds  to  Rusanov’s  scheme,  and  (4.5b")  is  the  Lax-Friedrichs 
scheme. 


The  following  simplification  allows  us  to  compute  a  single  “Riemann  solver"  per 
side,  rather  than  the  number  required  by  the  quadrature  formula  in  (4.3): 


(4.6a) 
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Here 


&k’*mL;R’dS 


’  “•  ~  |ac‘ 

r  /  f?,dS 

1  'ac- 

(4.6b) 

rfs  , 

~\ark\  f  fN(R*)dS  . 

\dcj\  ddc; 

(4.6c) 

Both  iij  and  /■  ran  be  computed  by  the  same  quadrature  formula  as  in  (4.3),  i.e., 


Uj  —  ^  ^  Qm^j  (^-m  i  V  ) 
m 

/‘  =  ^«m/4^mdn))  • 


(4.7a) 

(4.7b) 


We  observe  that  when  is  given  in  its  Taylor  series  form  (3.2),  it*  can  be  computed 
analytically  by 


r—  1 


*=0  |<|=fc 

=  jacpfi  Lc; (x  '  C;) 


(4.8a) 

(4.8b) 


which  is  usually  more  economical  than  the  computation  in  (4  7a).  The  expression  for 
i/(ui,u2)  in  (4.6a)  is  the  same  as  in  (4.5). 

Further  simplification  is  obtained  when  a  Taylor  expansion  of  /  (Rj(xjUn))  around 
Cj  is  available  to  us  in  the  form 


with 


Then,  as  in  (4.8), 


r_1  1 

f(R,(x-,u”))  =  Y,T,  £>-«)' fi  +  0(kr) 

*=o  *•  |/|=* 


«-^/(O  +  0(^M) 


/bEE  *f  (*>  • 

*  =  0  |f|=fc 


(4.9a) 


(4.9b) 


(4.10) 
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where  JV*  is  the  outward  normal  to  dCj.  Note  that  we  have  taken  only  a  single  value 
of  the  normal,  thus  for  boundaries  with  nonzero  curvature  the  RHS  of  (4.10)  has  to  be 
modified  accordingly. 

We  remark  that  in  the  constant  coefficient  case,  all  the  schemes  with  the  same 
v{ui,u-2)  are  identical  to  each  other,  and  the  scheme  with  v(u\,U2)  =  jA/v(ui, U2)|  is 
identical  to  the  original  scheme  (2.11). 


■  5.  Algorithmic  Considerations 


In  this  section  we  consider  the  simplest  scheme  of  the  previous  section 


dvj 

dt 

T 


-E 


f 

J  h* 


1  r 

2 


\ac; 


where 

=  E  E  b‘D’ 

k= 0  |£|  =  * 

=  EE  b>  (*>  ■ N!) 

Jt=0  |f|  =  Jt 

and  discuss  two  versions  which  differ  in  the  way  in  which  Ft  are  obtained. 


(5.1a) 


(5.1b) 


(5.1c) 


(o.  Id) 


In  the  first  version,  we  start  by  selecting  a  stencil  J(i)  to  the  cell  C,,  and  compute 
the  reconstruction  Rl(x;vn)  in  its  Taylor  series  form  (3.2).  Once  we  have  calculated 
{Dt}  we  evaluate  {Ft}  by  using  analytic  expressions  for 

F'=-^-'f(Ri(z;v"))Ucl  ■  (5.2) 

Still  within  the  cell  C,,  we  proceed  to  calculate  uf  and  /*  by  (5.1c), (5. Id)  respectively 
for  all  sides  k. 


After  doing  so  for  all  cells  Ct  in  P,  we  sweep  over  all  the  boundaries  dC,  to  compute 
the  numerical  flux  /,  „  by  (5.1b).  Finally,  we  go  over  all  the  cells  C{  to  evaluate  by 
(5.1a). 

This  version  makes  computational  sense  when  the  expressions  in  (5.2)  are  rela¬ 
tively  simple,  which  is  the  case  for  the  Euler  equations  of  compressible  perfect  gas  (see 
Appendix  2).  The  main  advantage  of  this  algorithm  is  that  most  of  the  computational 
work  is  done  within  the  cell,  with  minimal  communication  with  other  cells.  Hence  it 
seems  attractive  for  unstructured  grids,  especially  for  parallel  computers. 
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We  turn  now  to  describe  the  second  version  which  offers  a  low  operational  count 
at  the  cost  of  a  larger  storage  requirement  and  more  communication  between  cells.  The 
main  difference  from  the  previous  version  is  that  F(  are  computed  from  interpolation 
formula. 

We  start  this  algorithm  by  selecting  a  stencil  J(z)  for  the  cell  C{  and  compute  the 
reconstruction  Ri  (a-;  vn).  From  this  reconstruction  we  compute  the  point- values  at  the 
centroid  of  the  cell 

r— 1 

u?  =  Ri(c,-,v’)  =  v?-Y.'Ea‘.'D<  (5.3a) 

k=2  |/|=* 

and 

/”=/(«”)  •  (5.3b) 

Next  we  sweep  again  over  all  cells  to  compute 

=  (5-4a) 

Jh/(*;u")|,=Ci  .  (5.4b) 

Still  within  the  cell,  we  calculate  u*  and  /*  by  (5.1c)  and  (5. Id).  From  this  point 

on  we  proceed  exactly  as  in  the  previous  algorithm  to  compute  the  numerical  fluxes 
— k 

ft,*  by  (5.1b)  and  in  another  sweep  to  compute  the  RHS  of  (5.1a). 

The  notation  /(a;/n)  stands  for  the  point-value  approximation  (3.14),  and 
Jp-/ (a;  /n)  is  the  appropriate  coefficient  D( ,  which  is  computed  from  the  linear  sys¬ 
tem  (3.16).  We  point  out  that  the  same  set  of  cells  is  used  for  both  the  reconstruction 
step  and  the  interpolation  step  (see  the  discussion  at  the  end  of  Section  3);  this  is  of 
particular  importance  for  ENO  schemes  with  adaptive  selection  of  stencils.  Another 
observation  is  that  itj  can  be  computed  directly  from  the  reconstruction  using  D(  rather 
than  D(. 

This  second  scheme  has  a  low  operational  count,  especially  for  structured  grids 
where  the  computation  of  Df  and  Ft  is  accomplished  by  predetermined  finite  difference 
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expressions.  The  operational  count  for  this  scheme  is  1  selection  of  stencil  and  1  f(u) 
evaluation  per  cell,  and  1  “Riemann  solver’’  per  side.  We  shall  come  back  to  this  scheme 
in  Section  9  where  we  discuss  schemes  for  a  rectangular  grid. 

We  note  that  even  in  the  constant  coefficient  case,  the  two  versions  are  not  identical 
—  the  second  version  uses  many  more  cells  than  the  first  one. 


23 


6,  Fixed-Stencil  Schemes 


In  this  section  we  consider  the  schemes  described  in  Sections  4  and  5,  with  a  fixed 
central  stencil  J(i).  The  stencil  is  assigned  to  the  cell  on  the  basis  of  geometrical 
considerations  alone  so  that  it  will  be  as  centered  around  the  cell  as  possible.  Thus  in 
the  constant  coefficient  case 

f(u)  =  au  (6.1a) 

the  scheme  is  a  linear  operator  (unlike  the  ENO  schemes  of  the  next  section  where  the 
stencil  assigned  to  the  cell  depends  on  the  solution). 

We  note  that  although  the  stencil  for  the  reconstruction  step  is  centered,  the 
resulting  scheme  (2.7)  is  upwind  biased:  In  the  constant  coefficient  case 

E(T)R{x;vn)  =  R(x -aT-,vn)  (6.1b) 


and  therefore 

v”+i  =  A(Cj)R(x-aT-,vn)  .  ,  (6.1c) 

Based  on  analysis  of  some  simple  cases  and  some  numerical  experiments,  we  feel 
it  is  safe  to  conjecture  that  the  scheme  (2.7)  with  a  centered  stencil  J(i)  is  i^-stable 
in  the  constant  coefficient  case.  We  also  note  that  the  choice  of  centered  stencil  results 
in  the  most  accurate  reconstruction.  Therefore  such  schemes  are  excellent  numerical 
solvers  for  problems  with  smooth  solutions  or  with  weak  shocks;  in  the  latter  case  one 
can  add  some  form  of  high-order  numerical  dissipation  to  the  scheme,  e.g.,  one  can  add 
dissipation  in  the  form  of  a  filter  [12],  (10]. 

Clearly  when  the  stencil  J(i )  contains  a  discontinuity  of  the  solution,  the  recon¬ 
struction  Ri  has  spurious  oscillations  due  to  a  Gibbs-like  phenomenon.  In  this  case  we 
can  either  select  to  use  a  different  stencil  which  does  not  include  a  discontinuity  —  this 
technique  will  be  described  in  the  nexi  section,  or  to  modify  the  reconstruction  within 
the  same  stencil.  To  do  so  we  go  back  to  the  rather  old  ideas  of  hybrid  schemes  ([10], 
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[11])  except  that  here  we  hybridize  the  reconstructions  rather  than  numerical  fluxes. 
Symbolically  this  can  be  written  as 

Ri  ( x-vn )  =  6tv”  +  (1  -  6i)  R\  (x:vn)  (6.2a) 

where  r  >  2  is  the  high-order  reconstruction  (3.2),  and  v"  corresponds  to  the 
piecewise  constant  reconstruction 

Ri  (x;  vn)  =  wtn  ,  x  G  C,  (6.2b) 

which  is  monotone. 

The  automatic  switch  has  the  properties 

0  <  0,  <  1  ,  (6.3a) 

~  1  when  «7(i)  contains  a  discontinuity  ,  (6.3b) 

6,  —  0  (hr~1)  when  the  solution  is  smooth  in  J(i)  .  (6.3c) 

Rewriting  (6.2a)  as 

Rt(x]vn)  =  R1-  (x-,vn)  +  6i  [t>"  -  Ri(x;vn)]  (6.4a) 

we  see  that  since 

v? -Rrl(x-,vn)  =  0(h)  (6.4b) 

where  the  solution  is  smooth,  we  get  from  (6.3c)  that  the  formal  order  of  accuracy  is 
preserved. 

Preliminary  analysis  for  rectangular  grids  shows  that  one  can  construct  an  auto¬ 
matic  switch  which  satisfies  properties  (6.3),  and  that  the  resulting  scheme  is  essentially 
non-oscillatory.  We  refer  the  reader  to  Appendix  1,  where  we  present  analysis  for  the 
one-dimensional  case. 

In  this  paper  we  concentrate  on  deriving  ENO  schemes  by  an  adaptive  stencil 
technique  which  is  described  in  the  next  section. 
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7.  Adaptive-Stencil  Schemes 


In  this  section  we  consider  an  adaptive  reconstruction  technique  of  the  form  (3.2) 

r_1  1 

R(x;u  )  =  Ri  (x)u  )  =  ^2  X/  (x  ~  c,)(  De  ,  x  €  C,  (7.1a) 

*= o  '  |<|=|k 

D,  =  0  (c.)  +  0  (V-|<!)  (7.1b) 

r— 1 

■Dp  =  (cp, u  )  =  Uj  -  at>iDi  (7.1  c) 

fc=2  |<|=JL 

where  =  |^yJ4(C,)(x  -  c,)^ 

The  main  objective  in  this  adaptive  reconstruction  is  to  make  sure  that  if  C,  lies 
in  *he  smooth  part  of  the  function,  then  all  the  approximations  to  derivatives  D(  are 
computed  also  from  the  smooth  part  of  the  function.  This  guarantees  that 

ji2(x;  il)  -  u(x)l  =  0(hT)  for  x  6  C,  (7.2) 

for  all  cells  Ci  that  do  not  contain  a  discontinuity  and  thus  a  Gibbs-like  phenomenon  of 
spurious  oscillations  of  size  0(1)  in  the  neighborhood  of  a  discontinuity  is  avoided.  We 
remark  that  analysis  of  the  one-dimensional  case  in  [8]  shows  that  the  reconstruction 
is  generically  monotone  in  the  cell  that  contains  a  discontinuity. 

In  the  following  we  present  two  techniques  for  the  selection  of  appropriate  stencils, 
which  generalize  the  two  algorithms  that  have  been  used  in  the  one- dimensional  case 
in  the  earlier  development  of  the  ENO  schemes.  The  first  approach  is  to  consider 
several  candidate  stencils,  and  to  select  the  one  in  which  the  function  is  smoothest. 
The  second  approach  is  hierarchical:  We  begin  with  the  ith  cell,  and  at  each  step  of 
the  algorithm  we  add  another  cell  to  the  existing  stencil  for  the  computation  of  an 
additional  derivative. 

In  the  first  approach  we  have  to  specify  a  set  J  of  candidate  stencils.  Since  the 
choice  of  a  central  stencil  is  best  from  the  point  of  view  of  accuracy  and  stability,  let  us 
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agree  that  J  should  always  include  a  central  stencil  Jc(i),  which  should  be  selected  if 
the  function  to  be  reconstructed  is  smooth  in  it.  For  reasons  of  computational  efficiency 
we  would  like  the  number  of  stencils  in  J,  which  we  denote  by  |J|,  to  be  as  small  as 
possible.  On  the  other  hand  J  should  include  enough  directionally  biased  stencils  so 
that  we  can  select  a  stencil  in  which  the  function  is  smooth,  no  matter  where  we  place 
a  discontinuity  in  the  vicinity  of  Ct.  It  seems  to  us  that  for  each  side  dCk  of  the 
cell  we  need  to  have  a  directionally  biased  stencil  which  lies  in  a  conical  section,  the 
apex  of  which  is  the  cell  C,,  and  its  axis  of  symmetry  is  the  Inward  normal  to  dCk . 
This  stencil  widens  as  we  go  away  from  the  cell  Ci,  and  the  cone  is  truncated  once  it 
contains  enough  cells  to  determine  all  the  required  derivatives.  Let  us  denote  by  ,/*(?) 
the  diredtionally  biased  stencil  corresponding  to  the  side  dCk,  and  let  I\  denote  the 
number  of  sides  in  dCi,  so  the  minimal  set  of  candidate  stencils  seems  to  be 

J  =  ■  ■  ■ ,  •  (7.3) 

Thus  for  triangles  we  use  four  stencils  which  are  shown  schematically  in  Figure  1 . 


Let  us  denote  by  Dk  ,  |f|  >  1  the  approximation  to  derivatives  (7.1b)  which  is 
obtained  from  the  directionally  biased  stencils,  and  by  D\  the  values  that  correspond 
to  the  centered  stencil  Jc(i).  Let  ay 


Y,  ^'1  > 

ld=r-l 


(7.4a) 


serve  as  a  measure  of  smoothness  of  the  function  in  the  stencil.  Clearly  a  is  large  when 
the  stencil  contains  a  discontinuity;  if  the  function  is  smooth  in  the  stencil  then 


E 

Kl=r_ i 


<9V 

a?(Cl) 


(7.4b) 


Therefore  we  select  the  stencil  in  which  ay  is  minimal,  giving  preference  to  the  central 
stencil.  Following  Shu  [14]  this  can  be  done  by  using  aac  ,  |  <  o  <  1  for  the  central 
stencil  in  the  comparison. 


Next  we  would  like  to  describe  a  variant  of  the  above  technique,  which  selects  {Dt} 
in  (7.1b)  term  by  term,  without  selecting  a  single  stencil  for  the  reconstruction.  This 
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is  accomplished  by 


Dt  =  ma  (D\ ,  -  ■  ■  ,D(;Dcf)  ,  1  <  |£|  <  r  -  1  ,  -  <  a  <  1  (7.5a) 


Do  =  u.  -  ^ 

*=2  |f|  =  * 
r~1 

Rj(x-,u)  =  Dq  4-  (x  -  c,/  . 

fc=l  ’■  l<l=fc 


(7.5b) 


(7.5c) 


Here  ma  (xj ,  •  •  • ,  x*;  y)  is  a  modified  mininod  function  which  is  defined  by 


iQ(xi,--- ,xk;y)  =  | 


x,  if  m  =  |x,| 
y  if  m  =  a|yj 


(7.5a') 


where 


m  =  min  {  | x ,  | ,  •  •  • ,  |x*|,ct|j/|} 


We  turn  now  to  describe  the  hierarchical  algorithm  which  generalizes  the  one¬ 
dimensional  algorithm  in  [8]  to  multidimensions  and  general  geometries. 

We  denote  by  Jm(i )  the  stencil  of  m  cells  which  is  assigned  to  C,  at  the  mth  step 


of  this  algorithm 


Jm  (0  —  {*  1 1  '  '  ’  i  } 


(7.6a) 


and  by  J ^  the  indices  of  all  the  side- neighbors  of  Sm  =  (J  Dj,  i.e.,  the  cells  in  the 

j€.Jm 

exterior  of  Sm  which  share  a  common  side  with  dSm. 


With  Jm(i)  we  associate  an  invertible  system  of  m  linear  equations 


Cfm)d{rn)  =  u<m)  , 

rank  =  m 


(7.6b) 

(7.6c) 


which  is  obtained  from  writing  the  m  relations 


A(Cj)  R{{x\u)  =  Uj  ,  j  G  Jm{i) 
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in  the  form  (3.7) 


Qd  —  u 


(7.8a) 


where  Q  is  an  m  x  k (r)  matrix  and  d  =  (di,  •  •  •  ,  dK(r))T.  We  take  Q  *  to  be  the  first 
m  columns  of  Q,  d("d  is  the  first  m  components  of  d 


and 


d(m)  =  (d1,---,dmy 


u{m)  =  ,uim)T 


7.8b) 


(7.8c) 


We  begin  the  algorithm  by  setting 


(7.9a) 


then  for  m  =  1,  x(r)  —  1  we  define 

•^m-t-l(f)  =  Jm{i)  [^J  {fm+l}  i  *m+]  6  / 

In  order  to  select  we  consider  the  candidate  stencils 

J m  + 1  =  {*1  >  ’  "  '  ’  **»»  ij  }  '  J  ^ 


(7.9b) 


(7.10a) 


and  the  asociated  systems  of  (m  +  1)  linear  equations  of  the  form  (7.6b)  corresponding 
to  them 

_ (  «%X1  1  \  /  .  «  \ 

(7.10b) 


+  ,(m-f i)  (m  +  1)  •  j * 

Qj  d}  =  uj  ,  J  €  Jm- 


d{m+i) 

=  min 

j(m+D 

1m  + 1 

) 

Next  we  compute  d™+y  whenever  the  corresponding  system  is  invertible.  We  take  cm+j 
to  be  the  j  for  which  Id”14"1)  is  minimal,  i.e., 

(7.10c) 

i 

Here  |  |  denotes  some  weighted  norm  or  seminorm. 

After  completing  this  do  loop,  we  get  the  desired  reconstruction  (7.1)  by  setting 

d=d(K(r))  .  (7.10d) 
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We  observe  that  J ^  includes  all  the  side  neighbors  of  Jm(i)  and  thus  span  all 
possible  directions;  therefore  it  is  possible  to  find  j  in  J*n  so  that 


rank  +  =  m  +  1 


Furthermore,  if  Jm{i)  is  in  the  smooth  part  of  the  function  u(x)  and  J*t  includes  a  cell 
in  the  smooth  part  of  u(x )  which  is  directionally  suitable,  then  J„,  +  1(?)  will  also  be  in 
the  smooth  part. 

In  the  following  we  propose  an  efficient  implementation  of  this  algorithm,  which 
can  lie  viewed  as  an  ordered  Gauss  elimination  with  adaptive  row-pivoting. 

In  this  implementation  we  work  with  the  original  form  of  the  equations  (7.8a), 
which  are  ordered  as  follows:  In  the  mth  step  of  the  algorithm,  the  first  m  equations 
correspond  to  ?],•••,  im  in  Jm(i)  (7.6a).  These  are  followed  by  all  the  equations  corre¬ 
sponding  to  We  assume  that  at  the  beginning  of  the  step  the  first  m  equations 

are  in  upper  triangular  form,  and  that  di,  •  •  •  ,cfm-i  have  been  eliminated  from  the 
rest  of  the  equations.  As  is  customary  in  Gauss  elimination  we  add  the  RHS  of  the 
equations  as  an  extra  column  in  an  extended  matrix.  Thus  the  extended  matrix 
at  the  beginning  of  the  mth  step  is  as  follows 


Q(m)  = 


'(  T») 

" 

9n 

0 

~(m) 

Qmm 

0 

'■ 

(7.11) 


Jm- 1 


We  start  this  algorithm  by  setting  Q ^  to  be  the  1  x  (k  -j- 1)  matrix  which  corresponds 
to  the  equation  for  the  cell  C,. 

Given  we  show  now  how  to  evaluate  ?m+i  and  Q(m+1b  We  begin  by  adding 

to  Q(m)  the  equations  for  all  the  side-neighbors  of  C,m,  thus  completing  the  set  of  extra 
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equations  to  that  of  J ^n.  (Note  that  duplicity  of  equations  is  possible,  but  this  does 
not  interfere  with  the  execution  of  the  algorithm.)  We  use  the  diagonal  elements  in  the 
first  (m  —  1)  equations  to  eliminate  d ,,  •  •  •  .  dm_  i  from  the  additional  equations.  The 
next  stage  of  the  algorithm  is  to  use  qm'L  to  eliminate  dm  from  the  (m  +  l)th  equation 
and  on.  Doing  so,  we  are  now  in  a  position  to  form  Q1'™^ ,  which  is  an  upper  triangular 
form  of  the  system  (7.10b),  by  moving  any  of  the  equations  for  j  £  J ^  to  the  (m  +  1  )th 
row.  dK}  "r  ’  in  (7.10b)  can  now  be  computed  by  back  substitution  with  the  appropriate 
RHS  which  is  stored  in  the  (k  +  l)th  column  of  the  matrix. 

After  selecting  im +  i  by  (7.10c).  we  form  the  matrix  Q(m+1l  by  assigning  the  equa¬ 
tion  of  the  cell  zm+i  to  the  (771+  l)th  row.  This  completes  the  mth  step  of  the  algorithm. 

Once  we  have  computed  Q^K\  the  required  solution  d  is  obtained  by  back  substi¬ 
tution. 

We  remark  that  if  we  take  in  (7.10c) 

d(r',+1)  |  =  \dm+1\  (7.12) 

i.e.,  the  semi-norm  which  assigns  to  a  vector  the  absolute  value  of  its  last  component, 
then  there  is  no  need  to  back  substitute  in  every  step  of  the  algorithm.  Our  numerical 
experiments  seem  to  indicate  that  this  is  a  viable  practice. 

We  have  termed  this  procedure  as  ordered  Gauss  elimination  with  adaptive  row- 
pivoting  for  the  following  reasons:  "ordered"  -  because  we  feed  in  equations,  i.e.,  cells, 
in  a  particular  order;  “adaptive”  -  because  the  criterion  for  selecting  the  clustered  cells 
is  the  minimization  (7.10c)  of  the  derivatives  of  the  reconstruction  at  the  centroid  of 
the  ith  cell.  Within  the  context  of  adaptive-stencil  schemes  this  is  done  anew  at  each 
time  in  each  cell  and  the  result  depends  on  the  local  smoothness  of  the  reconstructed 
function. 

In  the  following  we  show  how  to  use  this  procedure  in  order  to  construct  fixed 
stencils  which  minimize  the  reconstruction  (or  interpolation)  error  of  a  given  function. 
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This  is  done  by  setting  the  problem  so  that  the  vector  d  in  (7.G)  and  (7.8)  is  the 
reconstruction  error  of  the  function  and  its  derivatives  at  the  centroid. 

Let  us  denote  by  f(x)  a  given  function  and  denote  its  derivatives  at  the  centroid 
of  the  ith  cell  by  f(. 


f  =  ~nci)  .  0  <  |<|  <  r  -  1 


and  its  cell-averages  by  f y 


f  j  =  A(Cj)  f 


(7.13a) 


7.13b) 


Let  us  define  the  vector  ( ,  r  =  e(f),  by 


k=0  |/|  =  jt 


(7.13c) 


and  consider  the  application  of  our  adaptive  algorithm  to  the  solution  of  the  set  of 


equations 


for  the  unknowns  Ef  =  E({f). 


Y,  ai’*E* = 

k  =  0  p|  =  Jt 


(7.14a) 


Let  us  denote  the  stencil  which  the  algorithm  assigns  to  C,  by  J(i\  /),  and  denote 
by  {F/}  the  solution  to  the  system 


H  a>'tFt  =  fj  '  J  €  ' 

k=0\(\  =  k 


(7.14b) 


clearly  {Fr}  is  an  approximation  to  {/^}  and 


E(  =  Ft  -  ft  . 


(7.14c) 


Let  us  denote  the  value  of  the  norm  in  (7.10c)  by  |F(/)|,  and  examine  the  way  in  which 
the  algorithm  arrives  at  the  stencil  J(i;  /):  We  start  with  the  cell  C,,  and  at  each  step 
of  the  algorithm  we  look  at  the  side  neighbors  of  the  existing  stencil  and  add  to  it  the 
one  which  minimizes  E(f).  We  observe  that  several  functions  / 1  (x),  •  •  • ,  fa(x)  can  be 
considered  simultaneously  by  adding  the  appropriate  terms  (7.14a) 
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as  an  extended  matrix  for  the  Gauss  elimination  and  use  in  (7.10c)  •  ,/s)| 

which  is  some  weighted  combination  of  ^(Z1  )[,•••.  |.E(/a)  | . 

This  observation  provides  a  very  useful  tool  for  the  construction  of  stencils  with 
special  properties.  Attempting  to  minimize  the  reconstruction  error  at  the  centroid  l. 
of  some  smooth  function  will  certainly  favor  a  centered  stencil;  if  we  use  a  discontinuous 
piecewise-smooth  function  f(x )  the  algorithm  will  select  a  stencil  which  is  as  centered 
as  possible  subject  to  the  constraint  that  it  should  not  include  the  discontinuity. 

The  simplest  way  to  construct  a  centered  stencil  of  n(r)  cells  is  to  apply  this 
procedure  to  monomials  of  degree  r 

f(x)  =  x *  ,  \e'\  -  r  .  (7.15a) 

We  note  that  e(f)  (7.13c)  takes  a  particularly  simple  form: 

t3  =  aJt(.  -  ]T  ah(  (7.15b) 

fc=0  \t\=k  Ci 

In  order  to  construct  the  directionally  biased  stencil  Jk(i )  for  (7.4)  in  an  automatic  way 
we  can  use  this  procedure  with  a  piecewise-polynomial  function  f{x)  of  degree  r  which 
is  discontinuous  at  the  face  dCf,  or  alternatively  to  use  a  combination  of  monomials 
of  degree  r  and  a  step-discontinuity  function  which  is  aligned  with  dCf . 

Another  way  to  construct  J*(t)  is  by  using  this  procedure  with  a  monomial  of 
degree  r,  but  restricting  the  side-neighbors  that  we  feed  into  the  algorithm  to  those 
which  are  contained  in  the  prescribed  conical  sections. 

We  remark  that  at  the  end  of  all  these  procedures  we  have  an  LU  decomposition 
which  is  needed  for  the  reconstruction;  as  is  customary  in  Gauss  elimination  we  store 
both  L  and  U  in  the  same  matrix. 

We  observe  that  this  technique  can  also  be  used  within  the  adaptive  stencil  schemes 
in  order  to  bias  the  selection  of  stencil  towards  a  central  one.  This  can  be  accomplished 
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by  adding  to  the  extended  matrix  (7.11)  a  column  of  e(/)  (7.13c)  for  some  smooth 
function  /,  and  properly  weigh  |d|  with  \E(f)\  in  (7.10c). 

We  remark  that  the  side-neighbors  which  we  feed  into  the  algorithm  are  restricted 
to  available  ones.  Thus  near  the  boundary  dV ,  the  growth  of  the  stencil  is  along  and 
away  from  dV. 

We  would  like  to  point  out  that  the  ENO  techniques  for  reconstruction  apply  as 
well  to  interpolation.  To  get  ENO  interpolation  all  we  need  is  to  replace  the  matrix  Q 
for  the  cell-averages  (3.7a)  by  the  matrix  Q  (3.16)  for  the  point  values. 

Finally  we  remark  that  up  to  now  we  have  stressed  the  desirability  of  biasing  the 
reconstruction  toward  a  central  stencil  for  reasons  of  accuracy  and  stability.  Taking 
into  account  the  cost  of  selecting  an  adaptive  stencil,  it  makes  sense  to  use  a  fixed 
central  stencil  altogether,  unless  it  contains  a  discontinuity.  In  order  to  decide  whether 
this  is  the  case,  we  can  use  the  quantity  8i  (6.3)  which  serves  as  an  automatic  switch 
in  the  context  of  hybrid  reconstruction  (6.4).  0,  is  0(1)  when  the  stencil  contains  a 
discontinuity  and  0(hr)  when  the  function  is  smooth  there.  Therefore  it  seems  possible 
to  determine  some  threshold  c,  so  that  an  adaptive  stencil  is  used  only  if  6i  >  c.  See 
Appendix  1  for  more  details. 
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8.  Hyperbolic  Systems  of  Conservation  Laws 


In  previous  sections  we  have  considered  the  reconstruction  of  a  piecewise-smooth 
function  n(x)  from  its  given  cell-averages  purely  as  a  problem  in  the  approximation 
of  functions,  without  actually  specifying  whether  u(x)  is  scalar  or  a  vector  function. 
When  we  are  dealing  with  a  vector  function  u(x)  on  the  level  of  approximation,  it 
is  up  to  us  whether  to  consider  its  components  as  independent  scalar  functions  or 
to  treat  the  whole  vector  as  a  single  entity.  Accordingly  we  can  either  apply  the 
adaptive  reconstruction  to  a  vector  function  in  a  component-wise  fashion,  which  means 
an  assignment  of  independent  stencils  to  each  of  the  components,  or  treat  the  vector 
as  a  single  entity  and  assign  a  single  stencil  for  all  components.  In  the  first  case  we 
use  |  |  in  the  expressions  of  Section  7  as  an  absolute  value  of  a  scalar  quantity;  in  the 

second  case  we  interpret  )  J  as  a  vector  norm. 

In  this  section  we  discuss  the  additional  aspects  that  come  from  the  fact  that  the 
reconstructed  function  u  is  a  solution  of  a  hyperbolic  system  of  conservation  laws. 

Since  u  —  u(x,t)  is  not  only  discontinuous  but  also  time-dependent,  it  is  possible 
for  discontinuities  to  come  close  to  each  other  and  interact.  Around  the  time  of  inter¬ 
action,  no  smooth  stencil  is  available  in  the  section  between  the  nearby  discontinuities 
and  some  spurious  oscillations  in  the  numerical  solution  can  be  generated.  Our  nu¬ 
merical  experiments  seem  to  indicate  that  the  component-wise  reconstruction  is  more 
robust  in  this  case  than  the  vector  reconstruction. 

Another  related  problem  is  the  fact  that  some  derived  quantities,  i.e.,  functions 
of  the  conserved  variable  u,  are  subject  to  constraints  imposed  by  the  physical  phe¬ 
nomenon  that  is  being  modelled,  e.g.,  density  and  pressure  are  nonnegative  quantities. 
Numerical  experiments  with  the  Euler  equations  of  compressible  gas  show  that  slight 
oscillations  in  the  conserved  variables  due  to  interaction  can  cause  much  larger  oscilla¬ 
tions  in  derived  quantities,  probably  due  to  the  fact  that  the  conserved  variables  are  out 
of  alignment.  Fortunately  we  find  that  this  problem  is  less  severe  in  multidimensional 
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calculations  than  it  is  in  the  one-dimensional  case,  where  all  interactions  are  necessarily 
frontal. 


In  the  one-dimensional  case  we  can  overcome  this  difficulty  by  using  locally  defined 
characteristic  variables,  since  this  set  of  variables  is  smoother  than  the  conserved  vari¬ 
ables  during  interaction  of  discontinuities.  This  technique  applies  in  a  straightforward 
manner  to  several  space  dimensions:  For  the  purpose  of  reconstructing  u  in  the  cell  C, 
we  consider  the  locally  defined  linear  transformation 

u=T(Tii)w  ,  io-T~1(v:)u  (8.1a) 

and  observe  that  due  to  the  linearity  of  the  transformation 


Wj  =  T  1  (ui)ilj  for  all  j 


(8.1b) 


We  apply  component- wise  reconstruction  to  the  values  in  (8.1b)  to  compute  at  the 
centroid  c,.  From  these  values  we  get  a  reconstruction  for  u  by  computing  derivatives 
of  u  at  the  centroid  by 

dlu  d(w 

8?  =  T(ui)a?e  '  (81c) 

C;  c, 


In  the  one-dimensional  case  we  take  T(u<)  to  be  the  matrix  of  eigenvectors  of  the 
Jacobian  |£.  To  get  a  smoother  set  of  locally  defined  characteristic  variables  in  several 
space  dimensions,  we  have  to  identify  the  normal  direction  to  the  discontinuity  in  u 
there,  and  use  the  eigenvectors  of  the  Jacobian  matrix  of  the  normal  flux.  We  remark 
that  since  (8.1)  is  completely  local,  the  formal  accuracy  of  the  reconstruction  does  not 
depend  on  the  choice  of  normal  direction. 


We  turn  now  to  describe  another  technique  that  works  well  for  solutions  of  Euler 
equations  of  compressible  gas  in  ID.  In  this  technique,  we  use  the  variables  ( p,q,p ),  i.e., 
density,  velocity,  and  pressure  in  order  to  reconstruct  the  conserved  variables,  which 
are  density,  momentum,  and  total  energy.  There  are  two  reasons  for  the  success  of  this 
approach: 
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(1)  p  and  p  are  variables  that  are  constrained  by  the  physics  of  the  problem  to  be 
nonnegative.  By  selecting  a  stencil  in  which  these  variables  are  smoothest  enables 
us  to  better  control  the  oscillations  in  these  quantities. 


(2)  q  and  p  are  continuous  across  a  contact  discontinuity,  which  is  the  central  wave 
in  the  Riemann  problem  corresponding  to  the  interaction.  Consequently  these 
variables  are  smoother  than  the  conserved  variables  in  regions  of  interaction. 

Let  us  now  return  to  the  general  case,  and  assume  that  there  is  a  preferred  set  of 
variables  iv(u)  which  is  a  nonlinear  transformation  with  a  well  defined  inverse  u(w). 
Let  us  also  assume  that  we  have  analytic  expressions  for  ^p-u>(u(x))  and  ~p-u(w(x)). 
The  quantities  immediately  available  to  us  are  w  (uj ),  which  by  (3.13b)  satisfy 

w(uj)  =  w(u(cj))  +  0  (h2)  .  (8.2a) 

The  coefficients  in  the  O  ( h 2)  term  above  involve  the  quantities 

AiCPix-c,)'  ,  |C  =  2  ,  (8.2b) 


which,  in  smoothly  varying  grids,  are  point-values  of  some  differentiable  function. 
Hence  for  grids  which  vary  smoothly  enough  we  can  use  a  component-wise  interpolation 
l(x\w(u))  in  (3.14)  to  get  0  ( h 2)  approximation  to  derivatives  of  w(n(x ))  at  the  cen¬ 
troid  c,.  From  these  values  we  can  compute 


c^u(u>(x)) 

dp 


1*1  =  1,2 


(8.2c) 


to  O  ( h 2)  and  thus  obtain  an  0  ( h 3)  reconstruction  Ri(x]u  )  which  results  in  a  third- 
order  accurate  scheme. 


We  observe  that  unlike  (8.1a)  the  quantities  {w (uj) }  are  defined  globally.  Hence 
this  procedure  for  third-order  schemes  is  easier  to  program  and  less  expensive  than 
that  of  the  locally  defined  characteristic  variables.  When  we  use  this  technique  for 
the  Euler  equations  of  gasdynamics  in  2D  and  3D  it  is  advisable  to  define  a  local 
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normal  direction  in  each  cell  and  represent  the  velocity  vector  in  normal  and  tangential 
components.  This  is  a  linear  transformation  which  does  not  interfere  with  the  formal 
order  of  accuracy  of  the  scheme. 


We  remark  that  this  technique  can  be  extended  to  the  general  case  in  two  ways: 


(1)  Use  the  values  of  jut*  (uj)}  for  the  purpose  of  selecting  a  stencil  in  which  w^(u ) 
is  smoothest.  Use  this  stencil  for  the  reconstruction  of  the  whole  vector  u,  and 
compute  derivatives  of  u  at  the  centroid  C{.  From  these  derivatives  of  u  compute 
analytically 

d(Wk(u) 

dxe 

*-  i 


to  the  desired  accuracy.  After  doing  so  for  all  components  A',  get  the  derivatives 


of  u  at  the  centroid  from  the  analytic  expression  for 


3u(u>(x) 


and  proceed  as  in  (8.1). 
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9.  Rectangular  Grids 


In  this  section  we  describe  the  schemes  of  previous  sections  for  the  case  of  rectan¬ 
gular  grids,  and  point  out  some  of  the  simplifications  that  result  from  considering  this 
case.  At  the  end  of  this  section  we  shall  describe  in  detail  the  third-order  case  which 
is  particularly  simple  and  seems  to  be  of  immediate  practical  importance.  To  simplify 
notations  we  consider  only  the  two-dimensional  case;  extension  to  3D  is  immediate.  As 
is  customary  we  denote  the  space  variables  by  x  and  y,  and  the  flux  components  by  / 
and  g. 

We  consider  the  IBVP 

ut+fx+gy  =  0  ,  ( x,y)eT>  ,  t>0  (9.1a) 

u(x,y,0)  =  u0(x,y)  ,  (x,y)el>  (9.1b) 

with  given  boundary  conditions  on  &D.  Typically  the  exterior  part  of  dV  is  rectangular. 
In  many  applications  V  contains  holes  which  correspond  to  rigid  objects,  in  which  case 
the  interior  part  of  dV  is  in  general  not  aligned  with  the  grid,  and  may  even  be  curved. 
In  the  latter  case  the  cells  which  are  side-neighbors  of  the  interior  part  of  dV  are 
treated  by  the  general  formulation  of  the  previous  sections.  In  here  we  shall  consider 
only  rectangular  ce'ls,  and  to  simplify  things  further  let  us  assume  that  the  grid  is 
uniform 


x,  =  ihx  ,  i  =  l,  •••,/*  ,  (9.2a) 

y3  =  jhy  ,  j  =  .  (9.2b) 

The  cells  are  identified  by  the  pair  (i,  j) 

C,j  =  [x„x,+  1}x[y},yJ+1}  ;  (9.3a) 

c,  j,  the  centroid  of  C,^,  is 

ci,j  ~  (^i+l/2iI/>  +  l/2)  •  (9.3b) 
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Discarding  the  multi-index  notation,  we  express  the  Taylor  expansion  of  u(x )  around 
Ci,j  by 


r— 1  ,  k 


(x,y)  =  h  12  ( * )  (*  -  x«+>/2)^  (v  ~  Vj+ i /2 ) 1 

k—0  K‘  t=0  '  ' 


■  0:r.(dyk~e 

Here  we  assume  that  u  is  sufficiently  smooth,  x  —  .r,+i/2  =  0  (hr).  y  —  Vj+i/2  =  0  (h  y). 
Applying  A(Cij)  to  u(x,y )  in  (9.4)  we  get  that  all  the  terms  with  odd  £  on  odd  k  —  £ 
in  the  summation  above  vanish  due  to  anti-symmetry;  thus  we  get 

Uij  =  a2k  [hid]  +  h2ydy)k  u|ft  .  +  0  (V^)  ,  (9.5a) 


—(  [x  «-f  1  /2 !  Vj+l/2  )  +  0{hr) 


where  [  ]  denotes  the  integer  part  and 


ak  =  2~k/(k  +  l)\  . 


(9.5b) 


For  r  =  6we  get  in  (9.5a) 


U{j  —  u  -f  24  {hxuxx  +  hyUyy) 

*b  jg20  zunxx  +  *-hxhyUXXyy  +  hyuyyyy )  ^  T  0(h  ) 


(9.5c) 


We  turn  now  to  describe  the  reconstruction  R  (x,  y;  u  ),  where  u  =  {ujj  }.  As  before 
we  rewrite  R{j  ( x,y,u  ),  the  polynomial  of  degree  (r  —  1)  for  the  cell  Cij,  in  the  form 
of  a  Taylor  expansion  around  the  centroid.  To  simplify  our  notation  we  translate  the 
origin  of  the  coordinate  system  to  the  centroid  and  scale  it  by  hx  and  hy  respectively, 


x  ~  xi+\j2 

hi 


y  -  yj+ 1/2 


for  simplicity  we  retain  the  notation  (x,y)  for  the  scaled  system. 


Rii(x,y,u)  =  D°'°  +  Y,T,Y,[c  )^V''D'‘-',  for  |x|  <  \  ,  |y|  <  \  .  (9.7a) 

*=1  K'  f=0  '  ' 
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The  accuracy  requirement  (3.5a)  can  be  expressed  by 


D{’k-(  =  (hXY(hy) 


k-i 


dku 


dx*dy 


k-e 


+  O  ( hr )  ;  0 < k  <r -  1  ,  0 < £ < k :  (9.7b) 


(0,0) 


the  terms  Dt,k  e  correspond  to  undivided  differences. 
In  an  analogous  way  to  (9.5),  we  get 

m  k 


k= i  f=o  v  ' 


,2  (k-i) 


(9.7c) 


which  is  equivalent  to  the  conservation  property  (3.1b).  For  r  =  6we  get  from  (9.7c) 
D°’°  =  ui}  -  ~  (D2’0  +  D0’2)  -  (D4’0  4  2D2-2  +  D0,4)  .  (9.7c') 

We  turn  now  to  describe  the  simple  form  that  the  numerical  scheme  (5.1)  takes 
for  rectangular  grids.  To  do  so  we  introduce  the  Taylor  expansion  of  the  flux  (4.9a)  in 
the  form 

=  F1'1  +  (5)  >V^'W  +  C(l’)  ,  (9.8a) 

*=rl  K"  t= 0  '  ' 


r— 1  „  k 
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=  G°'°  +  £  T,  £  (*)  x'yt~‘G>'t~‘  +  O (V)  (9.8b) 

*=1  K'  1=0  '  ' 


where  for  0  <  k  <  r  —  1  ,  0  <  £  <  k 

\£  (  l  \k—t 


=  (. h,y  (h,y 

Gtk-‘ =  (h,y  (h,)*-’ 
The  scheme  (5.1)  takes  the  form 


dk 


dxtdyk 
dk 


—t  f(u) 


dx(dyk 


~t  </(«) 


4-0(/ir) 

(0,0) 

+  0(hr) 


(9.8c) 

(9.8d) 


(0,0) 


dt  hx 

with  the  numerical  fluxes 


d  1  —  —  1 

—  U,)  ~  ~T~  {f t+l,j  —  / «,;)  —  7“  {9i,j+ 1  —  9t,>) 


(9.9a) 


fid  =  I  [/?-!,>  +  /£>  -  fi°>)  -  ".-id) 


(9.9b) 
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1 

9i'i  ~  2 

where  for  m  =  0, 1 


ffij- 1  +  9i,j  ~  t/(vij-i,Vij)  (Kj  ”  f;L-i) 


r"1  [k/2]  /  !-  i  1  \ 

=  ^°’° + £(-i)*(m+1w  E  ( 2* + 1 ) D 

*=1  f=0  '  ' 

r_1  /  >  I  1  \ 

=  -d0'0  +  E  u*  + 1 )  ° 

Jt  =  l  *=o  '  ' 

-1  lk/2]  /  t  \ 

fz  =  f°'° +^(-i)‘<’"+i)ot  £  ( *+ 1 )  f 

*=1  £=o  '  ' 

r-i  [*/2]  /  ,  ,  .  v 

iZ  =  G°'°  +  Y.  (  2e +  1  )  G 

k= 1  /=0  '  ' 


k-2t,2( 


2(,k-2( 


k- 2(,2( 


2(.k-2f 


ak  is  given  by  (9.5b). 


(9.9c) 


(9.9d) 


(9.9e) 


(9.9f) 


(9-9g) 


We  return  now  to  describe  an  algorithm  for  reconstruction,  i.e.,  how  to  compute 
Df  k~(  in  (9.7)  to  O  (hr)  from  the  given  cell  averages  {«,,>}•  The  most  convenient  way 
to  do  so  in  rectangular  grids  is  via  a  process  of  deconvolution  [9] ,  which  is  based  on  the 
observation  that  the  sliding-average  of  u 


u(x,y)  =  A  ( [a:  -  hx/2,x  +  hx/2]  x  [y  -  hy/2,y  +  hy/ 2])  •  u 
1  rh»/ 2  rh*/ 2 

=  TT~  /  u(x  +  +  Ti)d£dri 

nxny  J-ht/ 2  J-hx/2 


(9.10a) 


is  a  smooth  function  of  ( x,y )  (in  fact  it  is  one  order  smoother  than  u(x,y ))  and  that 
the  given  cell  averages  are  its  point  values  at  the  centroid,  i.e., 


Uij  =  u  (x,+1/2,y>+i/2)  (9.10b) 

Expanding  u(x  4-  £,  y  +  rj)  in  the  integral  as  a  Taylor  series,  we  get  that  the  relation 
(9.5)  holds  for  every  (x,y).  Therefore  we  can  differentiate  this  relation  to  express 
derivatives  of  the  sliding-average  u  in  terms  of  derivatives  of  u.  As  in  [9],  once  we 
properly  truncate  the  RHS  and  write  these  relations  as  a  system  of  linear  equations, 
we  get  an  upper  triangular  matrix  which  is  easily  inverted  by  back  substitution. 


43 


Before  deriving  this  system  of  equations,  we  introduce  the  notation 


r\ fc  — 

De,k~e  =  (Mf  (M*~*  a  t/k-t  +0(hr)  ,  0<fc<r-l  ,  0<(<A-  (9.11a) 

Ox  Oy  (Q  0, 


n°’°  - 

D  = 


(9.11b) 


We  derive  this  system  of  equations  by  symbolic  differentiation  of  (9.7c),  where 
differentiation  is  equivalent  to  increasing  the  appropriate  index.  Doing  so  we  get 


Dpq  P  =  Dp,q~p  + 


fr^t]  k 

it,  a2k  5Z  ( / )  D2t+pMk-t)+q-p , 

*=I  (=0  '  ' 


for  0  <  g  <  r  —  3  ,  0  <  p  <  r  —  3  —  q  , 

(9.12) 

We  remark  that  this  deconvolution  procedure  applies  also  to  smoothly  varying 
grids,  provided  that  the  corresponding  sliding-average  function  (9.10)  is  sufficiently 
smooth. 


Reconstruction  is  accomplished  by  substituting  an  appropriate  individual  differ¬ 
ence  for  jz)  ’  |  in  (9.12)  and  then  inverting  this  upper  triangular  system  of  linear 
equations  to  get  { Dm,n  }  by  back-substitution.  Thus  to  get  a  fixed-stencil  scheme  with 
a  centered  stencil  we  use  centered  undivided  differences  for  DP'9;  these  schemes  sire 
naturally  of  odd-order.  For  adaptive-stencil  schemes  we  use  differences  of  u  which  are 
computed  within  the  assigned  stencil.  We  observe  that  we  do  not  have  to  actually  se¬ 
lect  a  two-dimensional  stencil  and  that  the  same  computational  task  can  be  performed 
with  tensor  product  of  one-dimensional  stencils  as  follows: 

Step  1 

(i)  Using  values  of  m  =  0,  ±1,  ±2,  •  •  -,  apply  the  one-dimensional  algorithm 

(see  [8])  to  select  a  one-dimensional  stencil  of  r  cells  starting  with  i(i,j ),  and  use 
this  stencil  to  compute 

{d*;0}  ,  0  <  p  <  r  —  1  .  (9.13a) 
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(ii)  Similarly  in  the  y-direction,  select  a  y-stencil  of  r  cells  starting  with  j(i,j )  and 
compute 

{K'j9}  >  0  <  y  <  r  —  1  .  (9.13b) 

Step  2 


Compute  mixed  derivatives  DP'q  as  follows: 


(i)  Apply  Dpx ,  a  pth  order  finite- difference  operator  in  the  ^-direction  to  D 


n°'P 


DP'q  =  DID”"*  for 


rO,? 


r  —  1 


<  q  <  r  —  1  ,  1  <  p  <  min[y,  r  —  1  —  q]  (9.14a) 


using  values  of  D°’q  in  the  ar-stencil  for  (i,  j),  i.e., 


{A*!/}  >  0  <  i1  -i(i,j)  <  r  -  1 


(9.14b) 


(ii)  Similarly  in  the  y-direction  compute 


DP'9  =  DqyDP'Q  for 


<  p  <  r  —  1 


1  <  y  <  min(p,r  — 1  — p)  ,  (9.14c) 


using  values  of 

{-A\>°}  ’  0  <  J1  -K*J)  <  r  -  1  •  (9.14d) 

Mixed  derivatives  D *’*  that  are  evaluated  both  by  (9.14a)  and  (9.14b)  can  be 
either  averaged  or  minmod’ed.  The  values  used  in  (9.14b)  and  (9.14c)  should  be 
closest  to  (i,  j)  eis  possible. 


We  observe  that  due  to  averaging  over  cells  and  faces,  many  of  the  mixed  deriva¬ 
tives  cancel  out  in  the  expressions  for  the  point-value  and  the  numerical  flux.  As  we 
have  mentioned  in  Section  5,  the  quantities  {F*’*}  and  {<?*’*}  in  (9.9)  can  be  computed 
either  by  an  analytic  Taylor  expansion  (5.2)  using  values  of  {Dkt},  or  alternatively 
by  finite  difference  operators  applied  to  |/(Z?°’J°)|  and  |y(£)°)°)|.  We  can  compute 
these  derivatives  in  exactly  the  same  way  as  in  (9.13)-(9.14)  using  the  same  stencil. 
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We  turn  now  to  consider  the  case  r  —  3,  which  corresponds  to  a  particularly  simple 
scheme:  In  (9.12)  we  get 


00,0  = 


D°-2\ 


Dh0  =  DU°  ,  D2m  =  D2m 


D0,1  =  D0'1  ,  D°  2  =  D°’2 


in  (9.9)  vve  get  for  m  =  0. 1 


u™  =  D0  0  +  i(_l)m+1D,-°  4-  ^  (3D2,0  +  D0'2) 
=  F0,0  +  I(_i)m+  Iplfi  +  J_  ^3F2,0  +  F<>,2) 


(9.15a) 


(9.15b) 


(9.15c) 


(9.16a) 


(9.16b) 


and  symmetric  expressions  for  the  y-direction.  Note  that  we  do  not  have  to  compute  the 
mixed  derivative  term  D1'1  in  (9.15)  because  it  is  not  used  in  (9.16)  due  to  cancellation. 
Therefore  the  reconstruction  (9.12)  is  terminated  at  the  end  of  Step  1  (9.13)  for  the 
one-dimensional  derivatives.  Once  we  have  D2'0  and  D0’2  we  compute  D0,0  by  (9.15a), 


F°'°  =  /  ,  G°’°  =  g  (£>0,0)  . 


(9.17a) 


Now  jpi.o  F2,o  fo,2  50,1^0,2^2,0  can  ke  computed  either  by  an  analytic  expan¬ 
sion  (5.2)  or  by  differencing  the  values  of  F 0,0  and  G0,0  on  the  one- dimensional  stencils 
that  were  assigned  to  the  interpolation  of  u  in  Step  1,  i.e., 


ip2, 0  _  ip0,0  _  nrO.O  ,  rp0,0 

•’}  ~  i+2,j  «+l,j  ij  ’ 

F1  °  _  Fo,o  _  Fo,o  +  r  _  ■  _  2/2N  f2  0 
F0,2  =  F°’°  -  2F0,0  F°'° 

»-■>  .,j+ 2  .,>+1  + 

with  symmetric  expressions  for  G;  here  i  =  i(i,j)  and  j  =  j{i,j). 


(9.17b) 


In  Appendix  2  we  describe  the  implementation  of  this  third-order  accurate  scheme 
for  the  solution  of  the  Euler  equations  of  compressible  gas. 


46 


We  end  this  section  with  a  brief  review  of  the  fourth-order  scheme  (r  =  4).  From 
(9.12)  we  get 

P°'°  =u,-  ,  -  (d2'°  +D°'2)  ,  (9.18a) 


P1,  =  D 
P01  =  P 


o,i  _  n0>1  _  J_  Fn03  ,  n2’1 


(p°-3  +  p2>1)  , 


■»2 ,1  _  Tf2-1 


'\d0'2) 

1 

(9.18a) 

=  p2’0  , 

P3'0 

=  D30  , 

(9.18b) 

=  P°’2  . 

p°,3 

=  D"'3  , 

(9.18c) 

*  =  D1'3 

(9.18d) 

P1,0  +  — 
24 

(P3’° 

+  D1'2) 

,  (9.19a) 

in  (9.9)  we  get 

u™ ;  =  £>0,o  +  -1  (3D2’0  +  P0-2)  +  i(_l)"i+1  D1'0  +  ^(D3’0 +  D1’2)  ,  (9.19a) 

fm  =  Fo,o  +  J_  (3jF2,o  +  f0-2)  +  -(-1)",+1  [f]'°  +  —  (F3’0  +  F1’2)]  ,  (9.19b) 

24  v  2  24  ' 

and  symmetric  expressions  for  the  ^-direction.  After  completing  the  computation  of 

the  one-dimensional  derivatives  of  P  in  (9.13),  we  get  P 2,1  by  (9.14c),  i.e., 


°2.:!  =  ±  (o?:m.  -  -d?,;) 


(9.20a) 


depending  on  which  of  the  points  (i,j  ±  1)  is  included  in  the  y-stencils;  if  both  are  we 
2  1  • 

can  define  Di  '■  by  either 


Dij  =  min  mod  (f2’°  -  P2,°_i ,  P2;°+ 1  -  P2;°) 


(9.20b) 


D2’1  =  -  f D2,0  -  P2,0  ^ 

,  ]  2  V  ’  ■'+1  '-J"1  /  ’ 


(9.20b') 


12*  •  • 

is  computed  in  a  symmetric  way. 


When  we  compute  derivatives  of  /  and  g  by  differencing,  we  follow  the  exact  same 
procedure  that  was  used  for  computing  derivatives  of  u. 


47 


10.  Time  Integration 


The  semi-discrete  formulation  (4.1)  can  be  expressed  in  operator  form  as 


^=S(n(f))  ,  v{t)  =  {vj{t)} 

(10.1a) 

,  x  v-,  dC’f  -k 

S(V(t))J  =  Y,\Ck\ft,r  ■ 

(10.1b) 

In  order  to  construct  an  ENO  scheme  from  the  semidiscrete  formulation,  Shu  [1G]  has 

designed  a  class  of  multistep  time  integration  of  the  Runge-Kutta  (RK)  type  which  has 

the  property  that  when  applied  to  a  total- variation-bounded  (TVB) 

S,  it  is  also  TVB.  For  a  third-order  scheme  this  can  be  written  as 

spatial  operator 

v0  =  vn  ,  Kq  =  5  (v0) 

(10.2a) 

vj  =  i>o  -h  -f^o  ?  =  5(nj) 

(10.2b) 

i>2  =  +  -  (K0  +  K\)  ,  K2  =  5 (n2) 

(10.2c) 

=  o3  =  v"  +  l  (K0  +  Kx  +  4 K2)  . 

0 

(10. 2d) 

Since  vq,  Vi,  and  i>2  in  (10.2)  correspond  to  different  levels  of  time,  the  location  of 
discontinuities  in  these  functions  is  different.  A  discontinuity  in  the  solution  which  was 
located  near  the  boundary  of  the  cell  C*  for  no  may  very  well  be  located  at  a  neighboring 
cell  for  V\ .  This  observation  suggests  that  the  stencil- selection  procedure,  which  is  an 
expensive  part  of  the  algorithm,  has  to  be  repeated  at  each  step. 

In  the  following  we  present  a  single-step  time  integration  method 

’  (ia3a> 

—k 

which  is  a  modification  of  the  algorithm  in  [8].  The  numerical  flux  f  }  t  is  of  the  same 
form  as  (4.6a) 

7 h  =  .  (10.3b) 
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except  that  here 

fiK|ac‘|  l  LSEm‘]dsdt 

(10.3c) 

and 

“  t |ac*|  i  Lnm-W  - 

(10.3d) 

We  note  that  Uj(x,t)  =  E(t)Rj(x;vn)  is  the  small-time  solution  of  (2.1)  with  the 
smooth  polynomial  initial  value  u0(x)  =  R}(x;  vn ) .  The  approximations  that  we  use 
in  (10.3)  are  obtained  by  replacing  the  integrands  above  by  their  Taylor  expansion  in 
space  and  time  around  the  centroid  of  the  cell  c}  and  t  =  0: 


r-1  , 

Uj(x,t)  =  E{t)Rj  5Z  (X  ~  +0(hr) 


D 


k\ 

k= 0  |/|-fm  =  i 

d\(\  +  ' 


(<’m)  ~ 


dx\l  •  •  •  dxi‘  dtm 


u(x , t) 


x=Ci,t  =  o 


i — l 


J  ~  M  Z-,  V 

*=0  *"  P|  +  m  =  fc 

^K|  +  m 


r  =  Oi  ,t=0 


U, 


r— 1 


*=0  |<|+m=Jk 


r— 1 


fc=0  |^|  +  m  =  it 


'"Af.m)  +  0(hr)  , 

(10.4a) 

+  0(,,r-|<|-m)  . 

(10.4b) 

O(V)  , 

(10.5a) 

) 

(10.5b) 

form  as  (4.8a)  and  (4.10) 

^.m)  i 

(10.6a) 

nO'JV*)  ■ 

(10.6b) 

except  that  here 


6 


/t 

(/,m) 


(m  +  l)|€|!|dC* 


m  +  1 


(10.6c) 


and  bf  is  (4.8b);  t  =  (t\,  •  ■  •  ,£3 )  is  the  vector  of  s  indices  for  the  space  variables  and 
(e,m)  =  (*!,•••  ,  m)  is  the  extended  vector  of  m  -f  1  indices  where  the  last  one,  m, 

stands  for  time. 
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After  computing  {Z?(^0),  F((.o)  }  for  0  <  jfj  <  r  —  1  as  in  the  semidiscrete  scheme 
we  proceed  to  compute  D(i  m)  and  ,n)  for  1  <  m  <  r  —  1  ,  0  <  |^|  <  r  —  1  —  m  by  a 
Cauchy-Kowalewski  procedure  (see  [8])  as  follows: 


DO  for  m  =  1 ,  •  •  • ,  r  —  1 

Q\t\ +m  '  /  \ 

dx(dV,,U~  ^1V  ydr'dt”1-1  ^ ) 

^|f|  +  m  Qi  /  Qu  dmu\ 

d^dt^f  =  d^H  V"1 
END  DO. 


0  <  \i\  <  r  -  1  -  m 


0<|f|<r  —  1  —  m 


(10.7a) 

(10.7b) 


Here  we  have  used  differentiation  of  the  PDE  (2.1a)  to  get  (10.7a)  and  Hm  denotes 
the  functional  dependence 


dm 

dtrn 


f(u)  =  Hm 


(  du 


(10.7c) 


In  the  following  we  present  two  algorithms  for  the  implementation  of  (10.7)  which 
are  a  direct  extension  of  the  ones  described  in  Section  5  for  the  semidiscrete  case. 

Algorithm  1 

Compute  {D(<i0),F(<i0)}  ,  0  <  |f|  <  r  -  1  as  in  (5.2). 

DO  for  m  =  1,  •  •  ■ ,  r  —  1 

=  “  DIV  ,  0<|f|<r  —  1  —  m 

F((,m)  D,---,  Z?Km))  ,  0<K|<r-l-m  , 

END  DO. 


(10.7a') 

(10.7b') 


We  have  used  the  notation  “DIV”  for  the  analog  of  “div”  in  terms  of  indices;  /f/,m 
denotes  the  functional  dependence  in  (10.7b)  and  *  stands  for  spatial  derivatives  of 
order  less  or  equal  \t\. 
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We  observe  that  once  {D(/  0)}  are  computed,  all  other  calculations  in  (10.7')  are 
performed  within  the  cell. 

Algorithm  2 


Compute  {D(^o),F(;  o)}  ,  0  <  \t\  <  r  -  1  as  in  (5.4) 

DO  for  m  =  1 ,  •  •  • ,  r  —  1. 

(i)  For  all  j  in  the  computational  domain  compute 


£>(o,m)  =  -  DIV{.F(0)m_J)}  (10.8a) 

-F(0,m)  =  Hm  (D(0,0),  •  •  •  ,D(0lm))  (10.8b) 


(ii)  For  all  j  and  0  <  |^|  <  r  —  1  —  m  compute 


(10.8c) 


dl 

'(tm)  —  Qxf  I]\X'  F(0,m)J 


(10. 8d) 


END  DO. 


Here  J,(x;  •)  denotes  interpolation  on  the  stencil  that  is  assigned  to  the  cell  Cj. 

As  in  Section  5  we  observe  that  Algorithm  1  is  particularly  suitable  for  unstruc¬ 
tured  grids,  provided  that  the  analytic  expressions  for  Htm  are  available  and  are  sim¬ 
ple  enough;  this  is  the  case  for  the  Euler  equation  of  compressible  perfect  gas  —  see 
Appendix  2  The  second  algorithm  is  most  suitable  for  structured  grids  where  the 
differentiation  of  the  interpolation  in  (10. 8d)  can  be  expressed  by  some  differencing 
operator.  In  this  case  the  operational  count  for  Algorithm  2  is  rather  low. 

We  would  like  to  point  out  that  the  fully  discrete  schemes  described  above  differ 
from  the  abstract  scheme  (2.9)  in  one  important  feature:  They  use  E(t)  ■  Rj  (10.3d) 
rather  than  E(t)-R  in  (2.9).  The  use  of  the  solution  to  the  Riemann  problem  across  the 
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boundary  element  dCj  for  the  whole  time-step  ignores  the  interactions  at  the  vertices 
due  to  the  meeting  of  different  states  there.  Consequently  the  CFL  limitation  of  these 
schemes  is  rather  than  1  for  the  abstract  scheme  (2.9). 

Comparing  the  fully  discrete  versions  with  the  semidiscrete  formulation  it  seems 
clear  that  the  semidiscrete  formulation  is  easier  to  program  but  it  is  more  expensive  to 
use. 
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Appendix  1.  Hybrid  ENO  Reconstruction 

In  this  appendix  we  describe  some  preliminary  results  on  the  hybrid  reconstruction 
thai  was  symbolically  described  in  Section  6. 

First  some  notations:  We  consider  a  stencil  of  (r  +  1)  points  with  uniform  spacing 
h 

(u0,«i,...,Ur);  (Al.la) 

in  our  application  uj  are  taken  to  be  cell- averages.  Let  us  denote  by  T  the  translation 
operator 

THj  =  Uj+i  (Al.lb) 

and  by  A  the  undivided  forward-differencing  operator 


(Al.lc) 
and  denote 
(Al.ld) 


A  =  T  -  I,  I  =  T° 


A‘=(A)‘u,- 

Next  we  express  AJ  in  terms  of  {A*},  0  <  j  <  r  —  k  in  the  following  way 

a;  =  (Ar-*A{  =  (T  -  iy~k  A‘  =  E  fr  “ k)  (-i)'-1-'r'Aj 

(=0  '  ' 


Thus 


A;  =  EdAf,  a*  =  (-ir‘-'(r 

(=0  '  7 

We  define  6k ,  the  automatic  switch  for  the  fc-th  derivative,  by 


(A1.2) 


ek  = 


|A0r| 


t=  o 


r—  k 


r—k 


(A1.3) 


t=o 


t=o 


and  observe  that 


0  <  6k  <  1. 


(A1.4) 
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For  computational  puiposes  we  add  a  tolerance  e  to  the  denominator  in  ( A 1 . 3 )  and 
take  it  to  be  of  the  size  of  the  round-off  error. 


We  consider  now  two  important  cases 


I  A,*  |  =  6,,,„ 

A*  =  (-l)'lAfl, 


(Al.oa)(i) 

(A1.5b)(ii) 


where  6ftk  is  the  Kronicker-<5.  The  first  case  corresponds  to  a  step-discontinuity  in  u^k\ 
the  fc-th  derivative  of  u;  the  second  case  corresponds  to  a  discontinuity  in  u^m\  m  <  k. 
It  is  easy  to  see  that  in  both  cases  9k  =  1.  Thus 


6k  =  1  for  a  discontinuity  in  u^m\  k  >  m  >  0. 


(A1.5c) 


Next  we  consider  the  product  9kDk  in  the  case  that  u ^  is  continuous  and  as  in 
(9.11) 

Dk  =  hku(k)  +  0(hr)  (A1.6a) 


Clearly 


l«*5*l  =  O  (SfS)  •  [0(A*|u(i>|)  +  O(V)]  =  0(ftr);  (A1.8b) 


note  that  this  remains  so  even  when  it^  and  some  of  its  derivatives  vanish  at  a  point, 


u{k)  =  ...  =  u^+p-1)  =  o. 

In  this  case  the  denominator  in  9k  is  0(hk+p ),  but  so  is  Dk .  Thus 


(A1.7a) 


9k  =  0(hr~k~p),  9kDk=0(hr). 


(A1.7bj 


Next  we  consider  the  case  of  a  discontinuity  in  u^m\  k  <  m  <  r.  In  this  case 
AJ  =  0(hm)  and  A*  =  0(hk),  hence 


9k  =  0(hm~k),  Dk  =  hku(k)  +  0(hm) 


( A1.8 ) 
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We  summarize  all  these  cases  by 


(A1.9) 


{0  k  >  m  >  0 

Dk  +  0(hm)  k  <  m  <  r,  1  <  k  <  r  -  1; 
Dk  -f -  0(hr)  m  >  r 

here  m  is  the -order  of  first  discontinuous  derivative  of  it. 


We  turn  now  to  consider  the  r-th  order  reconstruction  (9.7)  in  the  one-dimensional 


=  D0+£pJJ*x‘,  |x|<i 

/c=l 

Dk  =  hku{k)(  0)  +  0(hr),  Dk  =  hku(k\ 0)  +  0{hr), 


(Al.lOa) 


(Al.lOb) 


where  { Dk }  are  obtained  from  {Dk  }  via  deconvolution,  i.e.  by  inverting  the  system  of 


linear  equations 


[X^l] 


£>«  =  £>*+  a2kD2k+\  0  <  g  <  r  —  1. 


( A 1 . 1 1 ) 


Setting  6°  =  0  we  define  the  hybrid  reconstruction  as  (A. 10)  in  which  {.D*}  are  obtained 
from  the  system 


‘  «  * 

(l-eq)Dq  =  Dq  +  «2 kD2k+q,  for  0  <  9  <  r  —  1. 


(A1.12) 


The  most  natural  choice  of  order  of  accuracy  for  this  hybrid  schemes  is  even  r  =  2  ? 
with  Dq  being  the  appropriate  central  differencing  for  the  stencil 


.,u0,...,u,),  r  =  2s. 


(A1.13) 


We  observe  that  if  u  has  a  discontinuous  derivative  in  the  stencil  (A. 13),  then 
from  (A1.9)  we  get  that  (1  —  9k)Dk  ■—  0  for  k  >  m  >  1  and  from  (A1.12)  Dk  =  0  for 
r  —  1  >  fc  >  m  >  1.  Hence  the  reconstruction  (Al.lOa)  becomes  a  polynomial  of  degree 


m  —  1. 
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We  turn  now  to  examine  the  case  r  =  2  of  a  second-order  scheme.  In  this  case 


D°  =  D°  =  «o, 
D1  =  (l-d')D1 


with 


Taking 


we  get 


i  =  lA-i  I  = 

|A$|  +  |Aj|  |A||  -f  |AI| 


JD1  =  i(u1-«_1)  =  i(A11  +  A‘) 


1[ia!i  +  ia;i-ia;-a;i] 


a|  +  a; 
lAii  +  iA"!' 


Recalling  that 


min(a,  b)  =  -  [a  +  b  -  |a  -  ft|] 


we  see  that 


(1  -01)/)1  =  /5‘min(lAoU^il)  if  sgn  (Aq)  =  sgn  (A})  =  5 
\  0  otherwise 

which  is  the  famous  minmod  function  m( AJ,  AJ).  Therefore 

R{x-,u)  =  u0  +  x  ■  m(Ao,  A}),  |x|  <  ^ 

which  is  TVD. 

For  r  =  4  we  get  in  (Al.12) 

D°  =  (1  -  e°)D°  =  D°  +  a2D2,  (1  -e^D1  =  P1  +  a2D\ 

(1  -62)D2=D2,  (1  -93)D3  =  D3 

which  is  inverted  to  give 

D3  =  (1  -63)D3,  D2=(1-62)D2,  D1  =  (1  -6')D'  -o2(l  -93)D3 
D°  =  D°  -a2(l  -02)D2 


(A1.14a) 

(A1.14b) 

(A1.14c) 

(A1.14d) 

(A1.14e) 

(A1.15) 

(A1.16a) 

(A1.16b) 
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with 


o'  =  |Ai2|/(|AL2K3|A1_1|  +  3|AJ|  +  |A!|). 

d2  =  |A12|/(|Ai2|  +  2|AiJ  +  |A2|),  (A  1.16c) 

o3  =  |A12|/(|A3.2|  +  iAij|), 

,  2  1 

D  =  “(U,  -  U- 1)  -  —  (t/2  -  U-j), 

4  j  5 

D2  =  -(uy  +  U-\)  —  +  U-2)  —  -Uq,  ( Al.lGd) 

1 

D  =  — (t/j  —  u_i )  +  ~(«2  —  U-2). 

The  extension  to  the  two-dimensional  reconstruction  (9.7)  is  straightforward: 


r—  1  k 


y-u)  =  d°'°  +  Y  ~  E  ( *  )  *fyk-(n('k-( 


k  =  1  f=0 


( Al.lTa) 


for  \x\  <  |y|  <  ^ 


with  Dp,q  given  by  inverting  the  system  of  linear  equations 


(1  -  0J)(1  -  8q-p)Dp’q~p  =  Dp’q-p+  Y  Q2*  Y  (  )  D2,+p’2U'-,)+9-? 

Jt=!  (=0^^ 

for  0  <  <7  <  r  —  1,  0<p<r  —  1  — 

(A1.17b) 

where  and  are  the  maximal  value  that  the  corresponding  one-dimensional  8k 
(A1.3)  takes  on  the  two-dimensional  stencil. 


Appendix  2.  Euler  Equations  of  Gas  Dynamics 


In  this  appendix  we  describe  a  particular  implementation  of  the  third-order  accu¬ 
rate  scheme  (9.15)-(9.1G)  to  the  Euler  equations  of  compressible  polytropic  gas: 


ut  +  fx  +  gy  —  0  (A2.1a) 

u  =  {p,pqT,pqy,E)T,  f(u)  =  qTu  +  (0,  P,0,  q1  P)1 ,  (A2.1b) 

g(u)  =  q»u  +  (0,0,P,q»P)T. 

with  the  equation  of  state 

P  =  (7-1)(E-  \pq\  (A2.1c) 

here  q2  =  ( qxf  +  (</y)2,  and  we  denote  mr  =  pqz.my  =  pqy . 

In  this  particular  implementation  we  use  the  primitive  variables 

w  =  w(u)  =  (p,  qz ,  qy,  P)T  ( A2.2) 


for  the  purpose  of  reconstruction  (see  the  discussion  in  section  8).  Given  cell-averages 
{uij}  we  compute  for  all  i,j 

w,i}  =  w(uij)  =  w(uij)  +  0(h2)  (A2.3a) 

and  apply  component-wise  selection  of  stencil  to  ia>,  i.e.  each  component  wk  is  assigned 
two  one-dimensional  stencils  ik(i,j),  jk(i,j).  Using  i  =  ik(i,j)  j  =  jk{i,])  we  compute 
for  1  <  k  <  4 


Dkt  =  (, hx)2wkxx\ij  =  w^2j  -  2u>f+1 ; 

Dk  =  hxwk\ij  =  wf+li.  -  Wk}  +  (*'-»'-  ^)£>rx 


and  similar  expressions  for  the  y-direction.  It  follows  from  (A2.3a)  that 


Dxz  =  (hx)2wrx(u)  +  0(hl),  Dz  =  hxwz(u)  +  0{h3x). 


(A2.3b) 

(A2.3c) 


( A2.3d) 
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(ii)  Calculate 

u  =  w  —  —  (i?IX  +  XJyy ),  /(u),  <7(u);  (A2.5a) 

note  that  u  is  accurate  to  0(h3),  and  therefore  also  f(u )  and  y(u). 

(iii)  Calculate  /x,  fxz,  fyy,gy,  gyy,gzx  from  the  analytic  expressions  of  derivatives 
of  f(u)  and  g(u),  i.e. 

/z  =  qzut  +  qzzu  +  (0  ,Px,Q,qzPx  +  qzzP)T, 

fxz  =  qzuZI  +  2qzzuz  +  +  (0,  P„  +  2(qz)xPz  +  qzZIP)T ,  (A2.5b) 

/yy  =  grUyy  +  2^'Uy  +  <7yyU  +  (0,  P yy,0,<?XP yy  +  2^X  Py  +  qyyP)T . 

Substituting  the  values  of  u  (A25a),  Dz,  Dzz,  Dy,  Dyy  (A2.3)  and  Dz,  Dxz,  Dy,  Dyy 
(A2.4)  in  the  RHS  of  (A2.5b)  we  get 

Fr  =  (hz)2f(u)z  +  O(hl),  FZI  =  (hz)2f(u)yy  +  0{h\ ),  (A2.5c) 
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FVy=(hl)f(u)„+0(hl), 
and  similar  expressions  in  the  y-direction. 

(iv)  Calculate  (9.16),  i.e.  for  m  =  0,1 


+  5(-l)'“+l^  +  ^(3Z>„  +  D„), 
/”  =  /(«.,)  +  +  Jj(3F„  +  F„), 


(A2.6) 


and  similarly  for  the  y-direction.  This  ends  the  calculation  within  the  cell. 

Next  we  compute  for  all  i.j  the  numerical  fluxes  by  (9.9b),  (9.9c).  After  doing 
that  we  can  form  the  RHS  of  the  semi-discrete  formulation  (9.9a)  and  update  one  of 
the  RK  steps  in  (10.2). 

As  we  have  mentioned  in  Section  10,  it  seems  rather  wasteful  to  repeat  the  cal¬ 
culations  in  (A2.3)-(A2.6)  and  (9.9)  three  times  for  all  the  steps  of  the  RK  algorithm 
(10.2),  when  we  can  complete  the  update  of  the  whole  time-step  by  using  (10.3)  with 
algorithm  1  in  (10.7')  i.e.  to  compute  a  modified  formula  for  um  and  /m  (10.6)  in 
which  we  add  time-derivatives  to  (A2.6) 

•*■!?  -  ■>,,  +  !(-ir+1(D,  +  (d,:)  +  \d,  +  lj(3 D„  +  D„)  +  \d„, 

1.  i  <A26') 


fZ  =  /(««)  +  5(-l)m+,(F,  +  jF„)  +  jF,  +  gj(3  F„  +  F„)  +  ±F„, 

and  similar  expressions  in  the  y-direction.  After  that  we  compute  numerical  fluxes  by 
the  same  formula  (10.6b)  using  a  single  Riemann  solver,  and  continue  to  compute  vn+1 
by  (10.6a).  To  do  so  we  have  to  modify  the  previous  algorithm  as  follows:  In  (A2.3) 
we  also  compute  mixed  space  derivatives  by 


=m(AiD^A^A^A*D*) 


(A2.3') 


where  m  is  the  minmod  function  and  A±  are  respectively  the  forward  and  backward 
undivided  difference  operators. 

In  (A2.4)  we  also  compute 


Djy  =  hZhyUXy  +  0^) 


(A2.4') 
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Ill  (A2.5)  we  also  compute 


FZy  =  hxhyf(u)xy  +  0(h 3),  Gzy  =  hzhy g(  u)xy  +  0{ h'5 ) 


(A?. 5') 


( A2.7a) 


After  (A2.5)  we  compute  time  derivatives  and  mixed  time-space  derivatives  by 

ut  =  -(/*  +  gy) 

uit  =  —{fxi  g j y) 
uyt  =  ~{fzy  +  9yy) 

ft  =qfa  +  q'ut  +  (0,Pt,0,qxP  +  qx  Pt)T 
fit  =  qztu  +  qx  ux  +  qxut  +  qxuxt 

+  (0>  Pxti  0,  qxtP  +  qtPx  +  qzPt  +  qT Pxt)T , 

and  similarly  for  gf,gyt; 

utt  =  —{fit  +  gyt ), 

jtt  =  qfu  +  2qf ut  +  q1  un  A  (0,  Ptt,0, qxtP  +  2qx Pt  +  qz Ptt)T- 
and  similarly  for  yit. 


(A2.7b) 


A2.7c) 


Note  that  we  do  not  compute  fxt  and  gyt  because  they  cancel  out  in  the  flux 
integration  and  thus  do  not  appear  in  (A2.6').  The  notation  that  we  use  in  (A2.6')  for 
time  and  space- time  derivatives  is 

Ft  =  rf(u)t  +  0(h 3),  Fxt  =  hxTf(u)xt  +  0(h 3),  Fu  =  t2  f{u)tt  +  0(h3)  (A2.7d) 
and  similarly  for  other  terms;  we  assume  r  =  0(h). 

We  remark  that  it  is  convenient  to  derive  the  quantities  Pt,  Pit,  Pyt,  Ptt  by  differ¬ 
entiating  the  equation  for  the  pressure 

Pt  +  qxPx  +  qyPy  +  7  P(qzx  +  qyy )  =  0,  (A2.8) 

and  derive  qf  ,qxt,qyt,qft  by  differentiation  of  the  relation 

mx  =  qxp, 

and  similarly  for  qy  (see  [8]). 

We  remark  that  for  purposes  of  reconstruction  in  (A2.3)  we  can  use  any  decom¬ 
position  of  the  velocity  vector  into  normal  and  tangential  components  rather  than  qz 
and  qy  (see  Section  8). 
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Appendix  3.  ENO  Schemes  in  Alternating  Dual  Grids 


In  this  paper  we  have  considered  a  set  of  cells  {Cj}  which  covers  the  computational 
domain  with  no  overlap.  We  have  introduced  a  reconstruction  procedure  R(x,u)  in 
which  we  assign  to  each  cell  a  polynomial  R}.  Time  evolution  is  done  by  solving  for 


small  time  the  IBVP  (2.8) 


wt  +  div  /  =  0 
w(x ,  0)  =  R(x]ii), 


(A3.1) 


followed  by  averaging  of  w(x,t)  on  the  same  set  of  cells. 


We  observe  that  for  a  time  step  which  is  limited  by  half  the  C'FL  number  the  value 
at  the  centroid  is  defined  in  terms  of  the  smooth  polynomial  problem 


wt  +  div  /  =  0 
w(x,Q)  =  Rj(x;  u). 


(A3. 2) 


Furthermore,  in  order  to  compute  the  cell-averages  we  have  already  prepared  the 
point  wise  space  and  time  derivatives  at  the  centroid,  so  the  evolution  of  the  centroid 
point- value  can  be  done  at  no  extra  computational  cost.  Hence,  we  can  easily  compute 
and  store  these  point-values  in  order  to  use  them  in  the  beginning  of  the  next  time 
level  for  purposes  of  reconstruction.  This  would  enable  us  to  get  pointwise  derivatives 
at  the  centroids  from  interpolation  of  these  point-values,  rather  than  from  direct  re¬ 
construction  of  averages  of  the  conserved  variables.  This  is  particularly  useful  if  we 
want  to  work  with  the  primitive  variables  for  purposes  of  selecting  smoother  data  for 
the  solution  of  Euler  equations  of  gasdynamics  as  discussed  in  Section  8. 

Once  this  is  done,  these  point-values  are  discarded;  the  point-values  are  computed 
anew  from  the  reconstruction  at  each  time  step.  Hence  it  is  a  side  calculation  for 
purposes  of  improved  reconstruction  and  the  point-values  themselves  do  not  constitute 
an  independent  set  of  variables. 

We  observe  that  the  point- values  at  the  centroid  do  not  enhance  the  reconstruction 
in  any  other  way  because  they  differ  from  the  cell-averages  only  by  an  0(h 2)  term  which 
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contains  second  and  higher  order  derivatives.  To  get  a  better  reconstruction  we  need 
the  two  sets,  point-values  and  averages,  to  be  separated  in  their  location.  Since  the 
only  reasonable  place  to  compute  point-values  is  in  the  centroid  of  the  cell  in  which  Rj 
is  defined,  we  look  for  the  possibility  of  averaging  in  a  different  kind  of  cell.  We  observe 
that  the  cell  formed  by  connecting  centroids  of  cells  around  a  vertex  is  a  suitable  one. 
We  refer  to  this  set  of  cells  as  vertex-cells  and  to  the  original  set  as  centroid  cells. 
When  we  deal  with  one  of  the  sets  we  refer  to  the  other  as  the  dual  one.  In  Figure 
A3  we  show  the  two  sets  of  cells  for  a  rectangular  grid  and  for  a  triangulated  mesh. 
We  observe  that  for  rectangular  grids  the  dual  sets  are  identical  except  for  a  shift;  in 
triangular  grids  the  two  sets  are  different:  the  centroid  set  is  made  of  triangles  while 
the  vertex  set  is  made  of  polygons  (typically  hexagons). 

We  use  these  two  sets  of  cells  in  an  alternating  fashion  in  time.  At  one  time-step  we 
assign  polynomials  Rj  to  one  set  of  cells  and  in  the  next  one  we  assign  them  to  the  dual 
set.  At  each  time-step  we  compute  point-values  at  the  set  at  which  the  polynomials 
Rj  are  defined,  i.e.  at  the  centroid  when  the  polynomials  are  defined  in  the  centroid 
set  and  at  the  vertex  for  the  time-step  in  which  the  polynomials  Rj  are  defined  in  the 
vertex  cells. 

We  turn  now  to  consider  the  computation  of  the  cell-averages.  Due  to  the 
divergence-free  form  of  the  PDE  the  averages  are  given  by  (2.8) 

\Cj\(v”+1  -  vj)  +  f  /  [/(«>(*,  0)  •  N}dS  =  0  (A3. 3) 

J0  JdCj 

where  w(x,t)  is  the  solution  to  (A3.1).  and  N  is  the  outward  normal  to  the  boundary 
dCr  When  we  consider  the  restriction  of  u;(a:,0)  to  the  cell  in  question  we  see  that 
the  middle  of  the  cell  (centroid  or  vertex)  is  an  apex  for  a  multiple  Riemann  problem. 
However,  each  side  of  the  cell  is  intersected  by  a  single  discontinuity.  Therefore  for  a 
time-step  which  is  limited  by  half  the  CFL  number,  the  computation  of 

/'  /  \{(w(x, t))-N\ds,  (A3. 4) 

Jo  J  dCj 

the  flux  across  the  side  dC involves  a  quasi-one  dimensional  Riemann  problem,  which 
is  formed  by  a  linear  segment  separating  two  smooth  functions  that  vary  in  space. 
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Fig.  A3.  Dual  cells  for  rectangular  and  triangular  grids 
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Since  these  are  smoothly  varying  functions,  the  solution  to  such  a  problem  can  be 
approximated  by  considering  a  linearized  one-dimensional  Riemann  problem  in  the 
direction  normal  to  the  linear  segment  (not  to  be  confused  with  the  normal  to  dCj );  the 
linearization  is  done  with  respect  to  the  jump  of  >r(.c.  0)  at  the  intersection  of  the  linear 
segment  with  dC *  using,  say,  Roe’s  technique.  We  demonstrate  this  approximation  in 
the  two-dimensional  rectangular  case.  First  let  us  align  the  coordinate  system  in  such 
a  way  that  the  discontinuity  is  along  the  x  axis  and  the  side  of  the  cell  is  along  the 
y- axis  between  y  =  —hj 2  to  y  =  hj 2.  We  consider  the  IVP 


wt  +  div  / 

«>(*,!/,  0)  = 


(  n'+(x. 


%y) 

%y)  y  <  o 


(A3. 5a) 


and  describe  how  to  approximate  the  numerical  flux  fT 

r  =  \  f  r 2  r(u'(0.y.t))dydt.  (A3. 5b) 

r'1  Jo  J-h/2 

To  do  so  we  use  Roe’s  linearization  with  respect  to  w  u  =  u>_(0,0)  and  =  te+(0,0), 
i.e.  we. consider  the  corresponding  constant  coefficient  case  with  a  matrix  A  defined  by 


fx(w°+)  -  fx(w°_)  =  -  w°_)  (A3. 5c) 


and  use  its  structure  to  approximate  the  integral  above.  Let  f±(y,t)  be  respectively 
an  approximation  to  the  flux  of  the  smooth  solution  of 


wt  +  div  /  =  0 
w(x,y,0)  =  u’±(x,y). 


(A3.6) 


We  approximate  f(w(0,y,t))  in  (A3. 5b)  by  fz(y,t) 


“f  <  V  <  M 

fZ(y,t)  =  \  fk(y,t)  akt  <  y  <  ak+it  k  =  \ . m-1 

.  U(y,t)  amt  <y  <  f 


where 


fk(y,t)  =  f-(y,t)  +  ^T[lj  ■  (/+  -  f-)]rj 


(A3. 7a) 


(A3. 7b) 


G7 


and  {aj,lj,rj}1jl=l  are  the  eigenvalues,  left-eigenvectors  and  right  eigenvectors  of  A. 
respectively.  Using  /z(y,t)  instead  of  fz(w(0,y,t))  in  the  integral  in  (A3. 5b)  we  get 
the  following  numerical  flux  for  this  side 


and 


/*  =  - 
J  rs 


</-+d)-£<w 
3=1 

2  r  /■“>* 


f_JU  -  f-)4ydt 


(A3. 8a) 
(A3. 8b) 


1  fT  fh/ 2 

f±  =  ~7  /  /  f±dydt.  (A3. 8b) 

r'i  io  J-h/2 


To  evaluate  these  integrals  we  use  the  Taylor  expansion  in  space  and  time  that  is 
described  in  Section  10.  We  remark  that  the  general  case  differs  from  the  above  only 
in  the  alignment  of  the  side  on  which  we  compute  the  numerical  flux,  i.e.  it  may  be 
skewed  and  uncentered. 


The  main  advantage  of  this  setting  is  that  we  double  the  number  of  values  that 
are  available  to  us  for  the  purpose  of  reconstruction;  this  effectively  doubles  the  spa¬ 
tial  resolution.  The  reconstruction  procedure  uses  now  a  combination  of  interpolating 
equations  (3.15)  and  averaging  equations  (3.6),  but  the  same  technique  of  adaptive 
selection  of  stencil  can  be  applied  to  this  combined  system. 

We  remark  again  that  the  set  of  point-values  is  not  an  independent  one  and  only 
plays  a  role  of  a  side  calculation.  The  method  can  be  viewed  as  a  scheme  for  cell- 
averages  in  which  we  alternate  the  set  of  cells.  The  most  natural  way  to  consider  the 
time  evolution  aspect  of  it  is  by  going  from  reconstruction  at  t„  to  reconstruction  at 

^n+l  • 

More  details  will  be  presented  in  a  future  paper. 


68 


NASA 

-Vo 

-H  <tvVj'v'<'Vi’,3IO 

1 .  Report  No 

NASA  CR-187637 
ICASE  Report  No.  91-76 
4  Title  and  Subtitle 


Report  Documentation  Page 


2. 


Government  Accession  No. 


3. 


5 


Recipient  s  Catalog  No. 


Report  Date 


MULTI-DIMENSIONAL  ENO  SCHEMES  FOR  GENERAL  GEOMETRIES 


September  1991 

6.  Performing  Organization  Code 


7 


.  Authorlsl 


8.  Performing  Organization  Report  No 


Ami  Harten 

Sukumar  R.  Chakravarthy 


9.  Performing  Organization  Name  and  Address 

Institute  for  Computer  Applications  in  Science 
and  Engineering 

Mail  Stop  132C,  NASA  Langley  Research  Center 
Hampton,  VA  23665-5225 _ 

12.  Sponsoring  Agency  Name  and  Address 

National  Aeronautics  and  Space  Administration 
Langley  Research  Center 
Hampton,  VA  23665-5225 


15.  Supplementary  Notes 
Langley  Technical  Monitor: 
Michael  F.  Card 


91-76 

10.  Work  Unit  No. 

505-90-52-01 

11.  Contract  or  Grant  No. 

NAS1-18605 

13.  Type  of  Report  and  Period  Covered 
Contractor  Report 

14  Sponsoring  Agency  Code 


To  be  submitted  to  Journal 
of  Computational  Physics 


Final  Report 
16  Abstract 


In  this  paper  we  present  a  class  of  ENO  schemes  for  the  numerical  solution 
of  multidimensional  hyperbolic  systems  of  conservation  laws  in  structured  and 
unstructured  grids.  This  is  a  class  of  shock-capturing  schemes  which  are  designed 
to  compute  cell-averages  to  high-order  accuracy.  The  ENO  scheme  is  composed  of  a 
piecewise-polynomial  reconstruction  of  the  solution  from  its  given  cell-averages, 
approximate  evolution  of  the  resulting  initial-value  problem,  and  averaging  of  this 
approximate  solution  over  each  cell.  The  reconstruction  algorithm  is  based  on  an 
adaptive  selection  of  stencil  for  each  cell  so  as  to  avoid  spurious  oscillations 
near  discontinuities  while  achieving  high  order  of  accuracy  away  from  them. 


17  Key  Words  (Suggested  by  Autborlsll 


18  Distribution  Statement 


ENO  schemes,  unstructured  grid, 
shock  capturing 


64  -  Numerical  Analysis 


Unclassified 

-  Unlimited 

i 

19  Security  Classif  (of  this  report! 

1  20  Security  Classif  (of  this  page! 

1  21  No  of  pages 

i  22  Price 

Unclassified 

j  Unclassified 

70 

:  A0  4 

NASA  FORM  1626  OCT  86 


NASA-Langley.  1991 


